"""
04-06  Normal Distribution Explorer
Chapters 04-02, 07-02

Run from the project folder:   python python/analysis.py
"""
import numpy as np
import pandas as pd
from scipy import stats
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

x = pd.read_csv("data/heights_cm.csv")["height_cm"].to_numpy()
n = x.size
mu, s = x.mean(), x.std(ddof=1)
N = stats.norm(mu, s)

# ── 1. TWO-WAY CALCULATOR ───────────────────────────────────────────────────
print("\n===== 1. NORMAL CALCULATOR =====")
print(f"Using the sample estimates: mean = {mu:.3f}, sd = {s:.3f}, n = {n}\n")

print("VALUE -> PROBABILITY")
for v in (160, 168, 175, 182, 190):
    print(f"  x = {v:3d}   z = {(v-mu)/s:+6.3f}   "
          f"P(X < x) = {N.cdf(v):.4f}   P(X > x) = {N.sf(v):.4f}")
print(f"\n  P(168 < X < 182) = {N.cdf(182) - N.cdf(168):.4f}")

print("\nPROBABILITY -> VALUE  (the inverse direction)")
for p in (0.01, 0.05, 0.25, 0.50, 0.75, 0.90, 0.95, 0.99):
    print(f"  P{100*p:<3.0f} = {N.ppf(p):7.2f} cm   "
          f"(critical z = {stats.norm.ppf(p):+6.4f})")
print(f"\n  middle 95%: {N.ppf(0.025):.2f} to {N.ppf(0.975):.2f} cm")

# ── 2. Z-SCORES ─────────────────────────────────────────────────────────────
print("\n===== 2. Z-SCORES =====")

z = stats.zscore(x, ddof=1)
print(f"CHECK  mean(z) = {z.mean():.10f}   sd(z) = {z.std(ddof=1):.10f}   "
      "(must be 0 and 1)")
print(f"|z| > 2 (unusual):      {(np.abs(z) > 2).sum()} of {n}  "
      f"({100*(np.abs(z) > 2).mean():.1f}%)")
print(f"|z| > 3 (very unusual): {(np.abs(z) > 3).sum()} of {n}  "
      f"({100*(np.abs(z) > 3).mean():.1f}%)")

top = np.argsort(-np.abs(z))[:5]
print("\nFive most unusual observations:")
print(pd.DataFrame({"height_cm": x[top], "z": z[top].round(3)}).to_string(index=False))

# ── 3. EMPIRICAL RULE CHECK ─────────────────────────────────────────────────
print("\n===== 3. EMPIRICAL RULE CHECK =====")

emp = pd.DataFrame({
    "k": [1, 2, 3],
    "lower": [round(mu - k * s, 2) for k in (1, 2, 3)],
    "upper": [round(mu + k * s, 2) for k in (1, 2, 3)],
    "observed": [round((np.abs(z) <= k).mean(), 4) for k in (1, 2, 3)],
    "empirical": [0.6827, 0.9545, 0.9973],
    "chebyshev": [np.nan, round(1 - 1/4, 4), round(1 - 1/9, 4)],
})
print(emp.to_string(index=False))
print("\nChebyshev is always satisfied but always weaker; the empirical rule is")
print("tight here because the data really is bell-shaped.")

# ── 4. NORMALITY ASSESSMENT ─────────────────────────────────────────────────
print("\n===== 4. NORMALITY ASSESSMENT =====")

print(f"mean {mu:.3f}   median {np.median(x):.3f}   difference {mu-np.median(x):.3f}")
print(f"skewness         {stats.skew(x):+.4f}   (0 if normal)")
print(f"excess kurtosis  {stats.kurtosis(x):+.4f}   (0 if normal)")
print(f"Pearson SK       {3*(mu-np.median(x))/s:+.4f}")

W, p_sw = stats.shapiro(x)
print(f"\nShapiro-Wilk  W = {W:.5f}, p = {p_sw:.4f}")
print("=> FAIL TO REJECT normality. The normal model is reasonable."
      if p_sw > 0.05 else
      "=> REJECT normality. Inspect the Q-Q plot before deciding what to do.")
print("\nRemember: with large n these tests reject trivial departures.")
print("The Q-Q PLOT is the primary evidence.")

# ── 5. NORMAL APPROXIMATION TO THE BINOMIAL ─────────────────────────────────
print("\n===== 5. NORMAL APPROXIMATION TO THE BINOMIAL =====")
print("n = 200, p = 0.55, find P(X >= 120)\n")

nb, pb = 200, 0.55
print(f"np = {nb*pb:.0f} >= 5  and  nq = {nb*(1-pb):.0f} >= 5   -> approximation is valid")
mu_b, sd_b = nb * pb, np.sqrt(nb * pb * (1 - pb))
print(f"mu = {mu_b:.1f}   sigma = {sd_b:.4f}\n")

exact = stats.binom.sf(119, nb, pb)
withcc = stats.norm(mu_b, sd_b).sf(119.5)
nocc = stats.norm(mu_b, sd_b).sf(120)

print(f"EXACT binomial              {exact:.6f}")
print(f"normal WITH correction      {withcc:.6f}   (error {abs(withcc-exact):.6f})")
print(f"normal WITHOUT correction   {nocc:.6f}   (error {abs(nocc-exact):.6f})")
print(f"\nThe correction reduces the error by a factor of "
      f"{abs(nocc-exact)/abs(withcc-exact):.0f}.")

# ── 6. PLOTS ────────────────────────────────────────────────────────────────
fig, axes = plt.subplots(2, 2, figsize=(11, 8.5))

axes[0, 0].hist(x, bins=20, density=True, color="#8A5FBF", edgecolor="white")
xs = np.linspace(x.min(), x.max(), 300)
axes[0, 0].plot(xs, N.pdf(xs), color="#0FA3A3", lw=2)
axes[0, 0].axvline(mu, color="#0B7A7A", lw=2, ls="--")
axes[0, 0].set(title="Heights with fitted normal curve", xlabel="Height (cm)")

stats.probplot(x, dist="norm", plot=axes[0, 1])
axes[0, 1].set_title("Normal Q-Q plot")

axes[1, 0].boxplot(x, vert=False, showmeans=True)
axes[1, 0].set(title="Boxplot", xlabel="Height (cm)", yticks=[])

xs = np.linspace(mu - 4*s, mu + 4*s, 400)
axes[1, 1].plot(xs, N.pdf(xs), color="#5B2A86", lw=2)
tail = xs[xs >= 185]
axes[1, 1].fill_between(tail, N.pdf(tail), color="#0FA3A3", alpha=0.5)
axes[1, 1].set(title="Shaded tail: P(X > 185)", ylabel="density")
axes[1, 1].text(mu + 2.4*s, N.pdf(mu) * 0.5, f"{N.sf(185):.4f}", color="#0B7A7A")

plt.tight_layout()
plt.savefig("normal_plots.png", dpi=110)
print("\nWrote normal_plots.png")