# Data setup: Load Titanic dataset for analysis
import warnings; warnings.filterwarnings('ignore')
import pandas as pd
import seaborn as sns
import numpy as np
np.random.seed(42)
df = sns.load_dataset('titanic')
print('Shape:', df.shape)
df.head()

# Probability: What fraction of passengers survived?
survived_prob = df['survived'].mean()
print('Probability a random passenger survived: {:.2f}'.format(survived_prob))

# Calculate basic statistics for age and fare
print('Age: mean={:.1f}, median={:.1f}, std={:.1f}'.format(df['age'].mean(), df['age'].median(), df['age'].std()))
print('Fare: mean={:.2f}, median={:.2f}, std={:.2f}'.format(df['fare'].mean(), df['fare'].median(), df['fare'].std()))

# Handling missing values safely
print('Missing ages:', df['age'].isnull().sum())
df = df.copy()
df['age_filled'] = df['age'].fillna(df['age'].median())
print('Nulls after fill:', df['age_filled'].isnull().sum())

# Simulate 10 coin flips, probability 0.5, 1000 times
import numpy as np
coin_flips = np.random.binomial(n=10, p=0.5, size=1000)
print('Simulated head counts (mean):', coin_flips.mean())

# Simulate number of siblings aboard to show Poisson distribution
from scipy.stats import poisson
sibs = df['sibsp']
print('Mean siblings:', sibs.mean())
pmf = poisson.pmf([0,1,2], mu=sibs.mean())
print('P(0): {:.2f}, P(1): {:.2f}, P(2): {:.2f}'.format(*pmf))

# Plot a histogram of age and compare to normal curve
import matplotlib.pyplot as plt
from scipy.stats import norm
ages = df['age_filled']
plt.hist(ages, bins=20, density=True, alpha=0.6, color='g')
xmin, xmax = ages.min(), ages.max()
x = np.linspace(xmin, xmax, 100)
p = norm.pdf(x, ages.mean(), ages.std())
plt.plot(x, p, 'k', linewidth=2)
plt.title('Age Distribution with Normal Curve')
plt.show()

# Exponential example: time between events can be modeled this way
from scipy.stats import expon
wait_times = expon.rvs(scale=2, size=1000)
plt.hist(wait_times, bins=30, color='c', alpha=0.7)
plt.title('Exponential Distribution of Simulated Wait Times')
plt.show()

# Simulate sampling fare 500 times, each sample size 30
sample_means = [df['fare'].dropna().sample(30, replace=True).mean() for _ in range(500)]
plt.hist(sample_means, bins=20, color='m', alpha=0.7)
plt.title('Sampling Distribution of Sample Means (Fare)')
plt.xlabel('Sample mean')
plt.show()

# Calculate 95 percent CI for the mean age
import scipy.stats as stats
ages = df['age_filled']
n = len(ages)
mean = ages.mean()
sem = ages.std(ddof=1) / n**0.5
ci = stats.t.interval(0.95, n-1, loc=mean, scale=sem)
print('95 percent CI for mean age: ({:.2f}, {:.2f})'.format(*ci))

# Test: Did women survive at a different rate than men?
survived_female = df[df['sex']=='female']['survived']
survived_male = df[df['sex']=='male']['survived']
t_stat, p_value = stats.ttest_ind(survived_female, survived_male)
print('T-statistic:', t_stat)
print('P-value:', p_value)

# Compute correlation and covariance for age and fare
corr = df['age_filled'].corr(df['fare'])
cov = df['age_filled'].cov(df['fare'])
print('Correlation:', round(corr,2))
print('Covariance:', round(cov,2))

# Fit and plot linear regression: fare ~ age
from sklearn.linear_model import LinearRegression
mask = df['fare'].notnull() & df['age_filled'].notnull()
X = df.loc[mask, ['age_filled']]
y = df.loc[mask, 'fare']
model = LinearRegression().fit(X, y)
plt.scatter(X, y, alpha=0.4)
plt.plot(X, model.predict(X), color='red')
plt.xlabel('Age')
plt.ylabel('Fare')
plt.title('Regression: Predicting Fare from Age')
plt.show()
print('Model slope:', model.coef_[0], 'Intercept:', model.intercept_)

# Logistic regression predicting survival from age and class
from sklearn.linear_model import LogisticRegression
df['pclass_filled'] = df['pclass'].fillna(df['pclass'].mode()[0])
X = df[['age_filled', 'pclass_filled']].dropna()
y = df.loc[X.index, 'survived']
clf = LogisticRegression(max_iter=200).fit(X, y)
score = clf.score(X, y)
print('Prediction accuracy:', round(score, 2))

# Try logistic prediction for a new passenger
input_age = float(input('Enter a passenger age: '))
input_pclass = int(input('Enter ticket class (1, 2, 3): '))
proba = clf.predict_proba([[input_age, input_pclass]])[0,1]
print('Predicted survival probability: {:.2f}'.format(proba))

# Load and plot airline passengers over time
url = 'https://raw.githubusercontent.com/jbrownlee/Datasets/master/airline-passengers.csv'
air = pd.read_csv(url)
air['Month'] = pd.to_datetime(air['Month'])
plt.plot(air['Month'], air['Passengers'])
plt.xlabel('Month')
plt.ylabel('Passenger Count')
plt.title('Monthly Airline Passengers')
plt.show()

# Estimate survival probability using Bayesian Beta
from scipy.stats import beta
success = df['survived'].sum()
fail = df['survived'].size - success
x = np.linspace(0, 1, 200)
prior_alpha = prior_beta = 1
posterior = beta(success + prior_alpha, fail + prior_beta)
plt.plot(x, posterior.pdf(x))
plt.xlabel('Survival Probability')
plt.ylabel('Density')
plt.title('Bayesian Posterior for Survival Probability')
plt.show()

# Simulate saving 200 people from the Titanic at random, 1000 times
pop_size = len(df)
boat_size = 200
num_trials = 1000
left_behind = [pop_size - np.random.choice(pop_size, boat_size, replace=False).size for _ in range(num_trials)]
plt.hist(left_behind, bins=15, color='orange')
plt.xlabel('People Left Behind')
plt.ylabel('Frequency')
plt.title('Simulation: Titanic Rescue')
plt.show()

# Group passengers into fare bands and calculate survival by band
df['fare_band'] = pd.qcut(df['fare'], 4, labels=['Low', 'Mid-low', 'Mid-high', 'High'])
grouped = df.groupby('fare_band')['survived'].mean()
print('Survival rate by fare band:')
print(grouped)