import warnings
warnings.filterwarnings("ignore")

# Import all libraries needed throughout
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from scipy import stats
np.random.seed(42)

# Simulate 10 coin flips (0 = tails, 1 = heads)
coin_flips = np.random.choice([0, 1], size=10)
print("Coin flips:", coin_flips)

# Calculate the proportion of heads
prop_heads = np.mean(coin_flips)
print("Proportion of heads:", prop_heads)

# Simulate 20 dice rolls (values 1 to 6)
dice_rolls = np.random.randint(1, 7, size=20)
print("Dice rolls:", dice_rolls)

# Probability of rolling a six
prob_six = np.mean(dice_rolls == 6)
print("Proportion of sixes:", prob_six)

# Example: describe dice rolls
print("Mean:", np.mean(dice_rolls))
print("Median:", np.median(dice_rolls))
print("Standard Deviation:", np.std(dice_rolls))

# Visualize the dice rolls
plt.hist(dice_rolls, bins=np.arange(1,8)-0.5, rwidth=0.8)
plt.xlabel('Roll')
plt.ylabel('Frequency')
plt.title('Histogram of 20 Dice Rolls')
plt.show()

# Introduce missing values on purpose
dice_with_nan = dice_rolls.astype(float)
dice_with_nan[::5] = np.nan  # Every fifth value becomes missing

# Check for missing data
print("Has missing?", pd.isnull(dice_with_nan).any())
# Fill missing values with the median
filled_dice = pd.Series(dice_with_nan).fillna(np.nanmedian(dice_with_nan))
print("Filled dice rolls:", filled_dice.values)

# Bernoulli: Simulate one coin flip (1 = heads, 0 = tails)
bernoulli_flips = np.random.binomial(1, 0.5, size=8)
print("Bernoulli results:", bernoulli_flips)

# Binomial: Simulate 10 trials of 5 coin flips each
binom_results = np.random.binomial(5, 0.5, size=10)
print("Binomial results (number of heads in 5 flips):", binom_results)

# Poisson: Simulate number of emails received per hour, avg 4, for 12 hours
poisson_counts = np.random.poisson(4, size=12)
print("Poisson counts (emails per hour):", poisson_counts)

# Draw histograms for each distribution
fig, axs = plt.subplots(1, 3, figsize=(13,4))
axs[0].hist(bernoulli_flips, bins=3, rwidth=0.8, color='skyblue')
axs[0].set_title('Bernoulli')
axs[1].hist(binom_results, bins=range(0,7), rwidth=0.8, color='orange')
axs[1].set_title('Binomial')
axs[2].hist(poisson_counts, bins=range(min(poisson_counts), max(poisson_counts)+2), rwidth=0.8, color='green')
axs[2].set_title('Poisson')
plt.tight_layout()
plt.show()

# Simulate 100 values from normal and exponential distributions
normal_data = np.random.normal(loc=0, scale=1, size=100)
exp_data = np.random.exponential(scale=1, size=100)

# Plot both distributions
fig, axs = plt.subplots(1,2,figsize=(10,4))
axs[0].hist(normal_data, bins=15, color='purple', alpha=0.6)
axs[0].set_title('Normal Distribution')
axs[1].hist(exp_data, bins=15, color='red', alpha=0.6)
axs[1].set_title('Exponential Distribution')
plt.tight_layout()
plt.show()

# Simulate population: waiting times (exponential, not normal)
population = np.random.exponential(2, 10000)

# Take 500 samples, each made of 30, and plot sample means
sample_means = [np.mean(np.random.choice(population, 30)) for _ in range(500)]
plt.hist(sample_means, bins=20, color='teal')
plt.title('Sample Means (CLT Demo)')
plt.xlabel('Sample Mean')
plt.ylabel('Frequency')
plt.show()

# Show actual population histogram for comparison
plt.hist(population, bins=20, color='grey', alpha=0.7)
plt.title('Original Population: Exponential Distribution')
plt.xlabel('Value')
plt.ylabel('Frequency')
plt.show()

# Data setup: Import the Titanic dataset using seaborn
titanic = sns.load_dataset("titanic")
print("Rows, columns:", titanic.shape)
print(titanic.head())

# Simulate sampling survival on the Titanic
samples = [titanic['survived'].dropna().sample(30, replace=True, random_state=42+i).mean() for i in range(300)]
plt.hist(samples, bins=20, color='navy')
plt.title('Sample Mean Survival Rates (n=30)')
plt.xlabel('Sample Survival Rate')
plt.ylabel('Frequency')
plt.show()

# Bootstrap confidence interval for Titanic survival rate
boot_means = [titanic['survived'].dropna().sample(200, replace=True).mean() for _ in range(1000)]
ci_low = np.percentile(boot_means, 2.5)
ci_high = np.percentile(boot_means, 97.5)
print(f"95% confidence interval: [{ci_low:.2f}, {ci_high:.2f}]")

# Permutation test for gender survival difference
women = titanic[titanic['sex']=='female']['survived'].dropna()
men = titanic[titanic['sex']=='male']['survived'].dropna()
actual_diff = women.mean() - men.mean()
combined = np.concatenate([women, men])
perms = [np.mean(np.random.permutation(combined)[:len(women)]) - np.mean(np.random.permutation(combined)[len(women):]) for _ in range(1000)]
p_value = np.mean(np.abs(perms) >= np.abs(actual_diff))
print(f"Actual difference: {actual_diff:.2f}")
print(f"Simulated p-value: {p_value:.3f}")

# Try with Titanic: Age and Fare
corr = titanic[['age', 'fare']].corr().iloc[0,1]
cov = titanic[['age', 'fare']].cov().iloc[0,1]
print(f"Correlation between age and fare: {corr:.2f}")
print(f"Covariance between age and fare: {cov:.2f}")

from sklearn.linear_model import LinearRegression
# Prepare the data
reg = LinearRegression()
X = titanic[['age']].fillna(titanic['age'].mean())
y = titanic['fare']
reg.fit(X, y)
print(f"Slope: {reg.coef_[0]:.2f}")
print(f"Intercept: {reg.intercept_:.2f}")

# Plot
plt.scatter(X, y, alpha=0.2, label='data')
plt.plot(X, reg.predict(X), color='red', label='fit')
plt.xlabel('Age')
plt.ylabel('Fare')
plt.title('Linear Regression: Fare vs. Age')
plt.legend()
plt.show()

from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
# Prepare features: age, sex, class
features = titanic[['age','pclass','sex']].copy()
features['age'].fillna(features['age'].mean(), inplace=True)
features['sex'] = features['sex'].map({'male':0, 'female':1})
labels = titanic['survived']
X_train, X_test, y_train, y_test = train_test_split(features, labels, test_size=0.2, random_state=42)
lr = LogisticRegression()
lr.fit(X_train, y_train)

accuracy = lr.score(X_test, y_test)
print(f"Test accuracy: {accuracy:.2f}")

# Data setup: Read air passengers data from URL
url = "https://raw.githubusercontent.com/jbrownlee/Datasets/master/airline-passengers.csv"
air = pd.read_csv(url)
print(air.head())

# Plot passengers by month
plt.plot(air['Month'], air['Passengers'])
plt.xlabel('Month')
plt.ylabel('Number of Passengers')
plt.title('Monthly Airline Passengers')
plt.xticks(rotation=45)
plt.tight_layout()
plt.show()

# Beta-Binomial example: coin with unknown chance of heads
prior_alpha, prior_beta = 2, 2  # Initial guess: any value is possible
data = [1, 0, 1, 1, 0, 1]  # Observed flips: 1=head, 0=tail
post_alpha = prior_alpha + sum(data)
post_beta = prior_beta + len(data) - sum(data)
x = np.linspace(0,1,100)
prior_pdf = stats.beta.pdf(x, prior_alpha, prior_beta)
post_pdf = stats.beta.pdf(x, post_alpha, post_beta)
plt.plot(x, prior_pdf, label='Prior')
plt.plot(x, post_pdf, label='Posterior', linestyle='--')
plt.xlabel('Probability of Heads')
plt.ylabel('Density')
plt.legend()
plt.title('Prior and Posterior Beliefs about Fairness of Coin')
plt.show()

# Simulate two coin tosses 1000 times
results = np.random.choice([0,1], size=(1000,2))
both_heads = np.sum(np.all(results==1, axis=1))
empirical_prob = both_heads / 1000
print(f"Number of both heads: {both_heads}")
print(f"Empirical probability: {empirical_prob:.2f}")

# Simulate survival rate CIs for Titanic classes
for cls in [1,2,3]:
    surv = titanic[titanic['pclass']==cls]['survived'].dropna()
    boot = [surv.sample(80, replace=True).mean() for _ in range(500)]
    lo, hi = np.percentile(boot, [2.5,97.5])
    print(f"Class {cls}: 95% CI [{lo:.2f}, {hi:.2f}]")
    

# Practice prompt: Try changing passenger class, or simulate class differences using passenger sex.
print("Practice: Repeat the above simulation but group by sex instead of class. What do you find?")