# Setup: Essential imports and settings
import warnings; warnings.filterwarnings("ignore")
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.tools.sm_exceptions import ConvergenceWarning, ValueWarning
warnings.filterwarnings("ignore", category=ValueWarning)
warnings.filterwarnings("ignore", category=ConvergenceWarning)
warnings.filterwarnings("ignore", category=RuntimeWarning)
 
 

# Coin flip basics: Run 10 flips
flips = np.random.choice(["H", "T"], size=10)
print("10 coin flips:", flips)
 
 

# Coin flips: Estimate the probability of heads
n_flips = 10000
results = np.random.choice(["H", "T"], size=n_flips)
prob_heads = np.sum(results == "H") / n_flips
print("Estimated probability of heads:", prob_heads)
 
 

# Dice rolling: Simulate 1000 rolls and show basic stats
rolls = np.random.randint(1, 7, size=1000)
mean_roll = np.mean(rolls)
std_roll = np.std(rolls)
print("Mean roll:", mean_roll)
print("Standard deviation:", std_roll)
 
 

# Missing data: Introduce and handle missing values
rolls_with_nan = rolls.copy().astype(float)
rolls_with_nan[::100] = np.nan
mean_nan = np.nanmean(rolls_with_nan)
print("Mean with missing data:", mean_nan)
 
 

# Simulate and plot a Binomial distribution
n_trials, p_successes = 10, 0.5
binom_sims = np.random.binomial(n=n_trials, p=p_successes, size=1000)
sns.histplot(binom_sims, bins=11, stat="probability")
plt.xlabel("Number of Heads in 10 Flips")
plt.title("Binomial Simulation (Coin Flips)")
plt.show()

# Simulate and plot a normal distribution
normal_data = np.random.normal(loc=0, scale=1, size=10000)
sns.histplot(normal_data, bins=40, kde=True, stat="density")
plt.title("Histogram of 10,000 Normal Samples")
plt.xlabel("Value")
plt.ylabel("Density")
plt.show()

# Central Limit Theorem: Sampling averages
pop = np.random.uniform(0, 10, size=100000)
sample_means = [np.mean(np.random.choice(pop, size=50)) for _ in range(1000)]
plt.hist(sample_means, bins=25, density=True)
plt.title("Sample Means (n=50) from Uniform Data")
plt.xlabel("Sample Mean")
plt.ylabel("Frequency")
plt.show()

# Confidence intervals: Estimate mean tip size
tips = sns.load_dataset("tips")
sample = tips["tip"].sample(n=30, random_state=42)
mean_tip = sample.mean()
std_tip = sample.std(ddof=1)
from scipy.stats import t
conf = 0.95
margin = t.ppf(1 - (1 - conf)/2, df=len(sample)-1) * std_tip / np.sqrt(len(sample))
lower, upper = mean_tip - margin, mean_tip + margin
print(f"95% confidence interval for mean tip: ({lower:.2f}, {upper:.2f})")

# Hypothesis testing: Are dinner tips bigger than lunch?
dinner_tips = tips[tips["time"] == "Dinner"]["tip"]
lunch_tips = tips[tips["time"] == "Lunch"]["tip"]
from scipy.stats import ttest_ind
stat, p = ttest_ind(dinner_tips, lunch_tips)
print("p-value:", p)
if p < 0.05:
    print("There is evidence that dinner tips are different from lunch tips.")
else:
    print("There is not enough evidence to say the tips are different.")
 
 

# Correlation and covariance: Are total bill and tip related?
corr = tips["total_bill"].corr(tips["tip"])
cov = tips["total_bill"].cov(tips["tip"])
print("Correlation:", corr)
print("Covariance:", cov)
 
 

# Simple regression: Predict tip from total bill
import statsmodels.api as sm
X = tips["total_bill"]
y = tips["tip"]
X = sm.add_constant(X)
model = sm.OLS(y, X).fit()
print(model.summary())

# Logistic regression: Predict Titanic survival
titanic = sns.load_dataset("titanic")
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
X = titanic[["pclass", "age", "sibsp", "parch", "fare"]].fillna(0)
y = titanic["survived"]
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=42)
logreg = LogisticRegression(max_iter=200)
logreg.fit(X_train, y_train)
score = logreg.score(X_test, y_test)
print(f"Accuracy of survival prediction: {score:.2f}")

# Time series basics: Airline passengers over time
url = "https://raw.githubusercontent.com/jbrownlee/Datasets/master/airline-passengers.csv"
airline = pd.read_csv(url, parse_dates=["Month"])
plt.plot(airline["Month"], airline["Passengers"])
plt.title("Monthly Airline Passengers (1949-1960)")
plt.xlabel("Date")
plt.ylabel("Number of Passengers")
plt.show()

# Bayesian intro: Update beliefs with Beta distribution
trials = 20
successes = 13
from scipy.stats import beta
x = np.linspace(0, 1, 100)
posterior = beta.pdf(x, successes + 1, trials - successes + 1)
plt.plot(x, posterior)
plt.title("Posterior after 13 successes in 20 trials")
plt.xlabel("Probability of Success")
plt.ylabel("Density")
plt.show()

# Mini-project Sim 1: Winning at dice
def play_game():
    rolls = np.random.randint(1, 7, 3)
    return np.sum(rolls) >= 15

results = [play_game() for _ in range(10000)]
prob_win = np.mean(results)
print(f"Probability of getting 15 or more in 3 rolls: {prob_win:.3f}")

# Mini-project Sim 2: Analyze real COVID-19 cases
import io
try:
    url = "https://raw.githubusercontent.com/owid/covid-19-data/master/public/data/jhu/new_cases.csv"
    covid = pd.read_csv(url)
    country = "United States"
    covid_long = covid.melt(id_vars="date", var_name="location", value_name="new_cases")
    cases = covid_long[covid_long["location"] == country]["new_cases"]
except Exception as e:
    # Fallback: Use random synthetic data if fetch fails
    np.random.seed(42)
    dates = pd.date_range("2021-01-01", periods=180)
    cases = pd.Series(np.random.poisson(20000, size=180), index=dates)
    cases.name = "new_cases"
    print("Warning: Could not load live COVID-19 data, using synthetic data for US instead.")
    cases = cases.reset_index(drop=True)
sim_samples = [np.mean(np.random.choice(cases.dropna(), size=14)) for _ in range(1000)]
plt.hist(sim_samples, bins=25, color="orange", alpha=0.7)
plt.title("Simulated Average COVID-19 Cases (2-week samples)")
plt.xlabel("Average New Cases per Day")
plt.ylabel("Frequency")
plt.show()