import pandas as pd
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

mtcars = pd.read_csv('mtcars.csv')
print(mtcars[['model', 'mpg']].head(3))

mpg = mtcars['mpg']
mu, sigma = stats.norm.fit(mpg)
print(f'Real mtcars mpg: n={len(mpg)}')
print(f'Fitted normal mean (mu): {mu:.3f}')
print(f'Fitted normal std dev (sigma): {sigma:.3f}')

x = np.linspace(mpg.min() - 2, mpg.max() + 2, 200)
pdf = stats.norm.pdf(x, mu, sigma)
plt.figure(figsize=(9, 4))
plt.hist(mpg, bins=10, density=True, color='steelblue', edgecolor='white', alpha=0.7, label='Real mpg histogram')
plt.plot(x, pdf, color='darkred', linewidth=2, label='Fitted normal PDF')
plt.title('Real mtcars MPG vs. a Fitted Normal Curve')
plt.xlabel('Miles per Gallon')
plt.legend()
plt.show()

within_1sd = mpg.between(mu - sigma, mu + sigma).mean()
within_2sd = mpg.between(mu - 2 * sigma, mu + 2 * sigma).mean()
print(f'Real share within 1 std dev of the mean: {within_1sd:.1%} (rule expects ~68%)')
print(f'Real share within 2 std devs of the mean: {within_2sd:.1%} (rule expects ~95%)')

aapl = pd.read_csv('aapl_ohlcv.csv').dropna(subset=['Close'])
daily_returns = aapl['Close'].pct_change(fill_method=None).dropna()
mu_r, sigma_r = stats.norm.fit(daily_returns)
print(f'Real trading days: {len(daily_returns)}')
print(f'Real mean daily return: {mu_r:.5f}')
print(f'Real std dev of daily return: {sigma_r:.5f}')

skew_r = stats.skew(daily_returns)
kurt_r = stats.kurtosis(daily_returns)
print(f'Real skewness: {skew_r:.3f} (0 = perfectly symmetric)')
print(f'Real excess kurtosis: {kurt_r:.3f} (0 = normal-like tails)')

plt.figure(figsize=(6, 6))
stats.probplot(daily_returns, dist='norm', plot=plt)
plt.title('Q-Q Plot: Real AAPL Daily Returns vs. Normal')
plt.show()

rng = np.random.default_rng(seed=11)
simulated_uniform = rng.uniform(low=0, high=10, size=5000)
print(f'Simulated sample mean: {simulated_uniform.mean():.3f} (theoretical: {(0 + 10) / 2:.3f})')
print(f'Simulated sample std: {simulated_uniform.std():.3f} (theoretical: {((10 - 0) ** 2 / 12) ** 0.5:.3f})')

plt.figure(figsize=(9, 4))
plt.hist(simulated_uniform, bins=30, density=True, color='steelblue', edgecolor='white', alpha=0.7, label='Simulated draws')
plt.axhline(1 / 10, color='darkred', linewidth=2, label='Theoretical uniform PDF (flat)')
plt.title('Simulated Uniform(0, 10): Draws vs. Theoretical Flat PDF')
plt.xlabel('Value')
plt.legend()
plt.show()

superstore = pd.read_csv('superstore_sales.csv')
superstore['Order Date'] = pd.to_datetime(superstore['Order Date'])
daily_counts = superstore.groupby(superstore['Order Date'].dt.date).size().sort_index()
busy_threshold = daily_counts.quantile(0.75)
busy_dates = pd.to_datetime(pd.Series(daily_counts[daily_counts > busy_threshold].index)).sort_values()
gaps_days = busy_dates.diff().dt.days.dropna()
print(f'Real busy-day threshold (75th percentile): {busy_threshold:.1f} order lines')
print(f'Real number of busy days: {len(busy_dates)}')
print(f'Real gap between busy days: mean={gaps_days.mean():.2f}, min={gaps_days.min():.0f}, max={gaps_days.max():.0f}')

loc, scale = stats.expon.fit(gaps_days, floc=0)
print(f'Fitted exponential scale (mean gap): {scale:.3f}')
x2 = np.linspace(0, gaps_days.max(), 200)
plt.figure(figsize=(9, 4))
plt.hist(gaps_days, bins=20, density=True, color='steelblue', edgecolor='white', alpha=0.7, label='Real busy-day gaps')
plt.plot(x2, stats.expon.pdf(x2, loc, scale), color='darkred', linewidth=2, label='Fitted exponential PDF')
plt.title('Real Gaps Between Busy Order Days vs. a Fitted Exponential')
plt.xlabel('Days Between Busy Days')
plt.legend()
plt.show()

p_within_2_days = stats.expon.cdf(2, loc, scale)
print(f'Modeled P(next busy day within 2 days) = {p_within_2_days:.4f}')
real_within_2 = (gaps_days <= 2).mean()
print(f'Real observed P(gap <= 2 days) = {real_within_2:.4f}')

fig, axes = plt.subplots(1, 2, figsize=(11, 4))
x3 = np.linspace(mpg.min() - 2, mpg.max() + 2, 200)
axes[0].plot(x3, stats.norm.pdf(x3, mu, sigma), color='steelblue')
axes[0].set_title('PDF: Fitted Normal for Real MPG')
axes[1].plot(x3, stats.norm.cdf(x3, mu, sigma), color='darkred')
axes[1].set_title('CDF: Fitted Normal for Real MPG')
plt.tight_layout()
plt.show()