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

rng = np.random.default_rng(22)

NAVY="#0a1230"; INK="#1a2138"; INK_SOFT="#4a5578"
CYAN="#0891b2"; PURPLE="#7c3aed"; AMBER="#d97706"; GREEN="#059669"; PINK="#db2777"; BLUE="#2563eb"; GRID="#e6e9f2"
plt.rcParams.update({"figure.facecolor":"white","axes.facecolor":"white","figure.dpi":110,"font.size":11,
   "axes.edgecolor":GRID,"axes.grid":True,"grid.color":GRID,"axes.axisbelow":True,"axes.spines.top":False,
   "axes.spines.right":False,"axes.titlesize":12.5,"axes.titleweight":"bold","axes.titlecolor":INK,"legend.frameon":False})
print("✅ Environment ready.")

✅ Environment ready.

income = rng.lognormal(mean=10.5, sigma=0.7, size=4000)   # right-skewed, all positive
logged = np.log(income)

fig, (a1, a2) = plt.subplots(1, 2, figsize=(12, 4))
a1.hist(income/1000, bins=60, color=AMBER, edgecolor="white")
a1.set_title(f"Raw income   (skew = {stats.skew(income):+.2f})"); a1.set_xlabel("$ thousands"); a1.set_yticks([])
a2.hist(logged, bins=60, color=GREEN, edgecolor="white")
a2.set_title(f"log(income)   (skew = {stats.skew(logged):+.2f})"); a2.set_xlabel("log income"); a2.set_yticks([])
plt.tight_layout(); plt.show()
print("Skew fell from strongly positive to near 0. On the log scale, equal distances mean equal RATIOS (percent change).")
print("Zeros or negatives? log breaks. Use log1p, a shift, or Yeo-Johnson (Demo 4).")

Skew fell from strongly positive to near 0. On the log scale, equal distances mean equal RATIOS (percent change).
Zeros or negatives? log breaks. Use log1p, a shift, or Yeo-Johnson (Demo 4).

x = rng.lognormal(mean=1.0, sigma=0.6, size=3000)   # right-skewed, positive
rungs = [("x^2 (up)", x**2), ("x (none)", x), ("sqrt(x) (down)", np.sqrt(x)), ("log(x) (down)", np.log(x))]

fig, axes = plt.subplots(2, 2, figsize=(11, 7))
for ax, (name, vals) in zip(axes.ravel(), rungs):
    ax.hist(vals, bins=40, color=CYAN, edgecolor="white")
    ax.set_title(f"{name}    skew = {stats.skew(vals):+.2f}"); ax.set_yticks([])
plt.tight_layout(); plt.show()
print("Up the ladder (x^2) makes the right skew WORSE; down the ladder (sqrt, then log) pulls it toward symmetry.")

Up the ladder (x^2) makes the right skew WORSE; down the ladder (sqrt, then log) pulls it toward symmetry.

low  = rng.poisson(4,  2000)     # small counts
high = rng.poisson(80, 2000)     # large counts
print("RAW counts:")
print(f"  low  group: mean {low.mean():5.1f}, variance {low.var():6.1f}")
print(f"  high group: mean {high.mean():5.1f}, variance {high.var():6.1f}   <- variance tracks the mean")
print("\nAfter sqrt:")
print(f"  low  group: variance {np.sqrt(low).var():.2f}")
print(f"  high group: variance {np.sqrt(high).var():.2f}   <- now both near ~0.25, stabilized")

fig, (a1, a2) = plt.subplots(1, 2, figsize=(12, 3.8))
a1.hist(low, bins=20, alpha=0.7, color=CYAN, label="low"); a1.hist(high, bins=30, alpha=0.7, color=AMBER, label="high")
a1.set_title("Raw counts: spreads differ wildly"); a1.set_yticks([]); a1.legend()
a2.hist(np.sqrt(low), bins=20, alpha=0.7, color=CYAN, label="low"); a2.hist(np.sqrt(high), bins=30, alpha=0.7, color=AMBER, label="high")
a2.set_title("After sqrt: comparable spreads"); a2.set_yticks([]); a2.legend()
plt.tight_layout(); plt.show()

RAW counts:
  low  group: mean   4.0, variance    4.2
  high group: mean  80.0, variance   81.4   <- variance tracks the mean

After sqrt:
  low  group: variance 0.32
  high group: variance 0.25   <- now both near ~0.25, stabilized

x = rng.lognormal(mean=2.0, sigma=0.8, size=3000)
bc, lam = stats.boxcox(x)        # returns transformed data and the ML-optimal lambda
print(f"Box-Cox chose lambda = {lam:.3f}  (near 0 -> essentially a log)")
print(f"skew: raw {stats.skew(x):+.2f}  ->  Box-Cox {stats.skew(bc):+.2f}")

fig, axes = plt.subplots(1, 3, figsize=(14, 3.8))
axes[0].hist(x, bins=50, color=AMBER, edgecolor="white"); axes[0].set_title("Raw (right-skewed)"); axes[0].set_yticks([])
axes[1].hist(bc, bins=50, color=GREEN, edgecolor="white"); axes[1].set_title(f"Box-Cox (lambda={lam:.2f})"); axes[1].set_yticks([])
stats.probplot(bc, plot=axes[2]); axes[2].set_title("Q-Q plot: close to the line = normal")
axes[2].get_lines()[0].set_color(GREEN); axes[2].get_lines()[1].set_color(PINK)
plt.tight_layout(); plt.show()
print("In ML, fit the transform (the lambda) on TRAIN only, then apply to test, same no-leakage rule as scaling.")

Box-Cox chose lambda = -0.000  (near 0 -> essentially a log)
skew: raw +3.49  ->  Box-Cox -0.00

In ML, fit the transform (the lambda) on TRAIN only, then apply to test, same no-leakage rule as scaling.

x = rng.lognormal(mean=3.0, sigma=0.9, size=5000)
arith = x.mean()
geo   = np.exp(np.log(x).mean())     # back-transformed mean of logs = geometric mean
med   = np.median(x)
print(f"arithmetic mean      : {arith:8.1f}")
print(f"exp(mean of logs)    : {geo:8.1f}   <- the GEOMETRIC mean, biased LOW")
print(f"median               : {med:8.1f}   <- back-transforms cleanly (monotone preserves order)")

fig, ax = plt.subplots(figsize=(9, 4))
ax.hist(x, bins=70, color="#c7ccda", edgecolor="white")
ax.axvline(arith, color=PINK,  lw=2.5, ls="--", label=f"arithmetic mean {arith:.0f}")
ax.axvline(geo,   color=GREEN, lw=2.5, ls="-",  label=f"exp(mean log) = geometric {geo:.0f}")
ax.axvline(med,   color=BLUE,  lw=2.5, ls=":",  label=f"median {med:.0f}")
ax.set_xlim(0, np.percentile(x, 99)); ax.set_yticks([]); ax.set_xlabel("value"); ax.legend()
ax.set_title("Back-transforming the mean of logs lands below the arithmetic mean")
plt.tight_layout(); plt.show()

arithmetic mean      :     29.8
exp(mean of logs)    :     20.1   <- the GEOMETRIC mean, biased LOW
median               :     20.0   <- back-transforms cleanly (monotone preserves order)

# --- Real-World beat: compare transforms by how far they cut the skew ---
from scipy import stats
BASE = "https://raw.githubusercontent.com/johnfisher-ai/Statistics-Data-Science-AI-Visual-Book/main/data/"
try:    api = pd.read_excel("../../data/transformations--response_times.xlsx", sheet_name="Data")
except FileNotFoundError: api = pd.read_excel(BASE+"transformations--response_times.xlsx", sheet_name="Data")
r = api.response_ms
bc, lam = stats.boxcox(r)
variants = {"raw": r, "log": np.log(r), "sqrt": np.sqrt(r), f"Box-Cox (lam={lam:.2f})": pd.Series(bc)}
for name, v in variants.items():
    print(f"{name:22s} skewness = {pd.Series(v).skew():+.2f}")
print("\nlog and Box-Cox pull the skew to about 0; the square root only helps partway.")
fig, ax = plt.subplots(1, 4, figsize=(14,3.4))
cols = [CYAN, GREEN, AMBER, PURPLE]
for a, (name, v), c in zip(ax, variants.items(), cols):
    a.hist(v, bins=30, color=c, alpha=0.85, edgecolor="white")
    a.set_title(f"{name}\nskew {pd.Series(v).skew():+.2f}", fontsize=10); a.set_yticks([])
plt.tight_layout(); plt.show()

raw                    skewness = +1.66
log                    skewness = +0.10
sqrt                   skewness = +0.84
Box-Cox (lam=-0.06)    skewness = +0.00

log and Box-Cox pull the skew to about 0; the square root only helps partway.

#	Demo	Idea it builds
1	Log tames the tail	a right-skewed variable to near-symmetric
2	The ladder of powers	one variable, four rungs, four shapes
3	Square-root for counts	stabilize variance on Poisson data
4	Box-Cox picks the power	maximum-likelihood lambda + Q-Q check
5	The back-transform trap	exp(mean of logs) is the GEOMETRIC mean

🎯 What you'll build in this notebook¶

⚙️ Setup, imports & the book's plotting style¶

⏱️ Real-World Example: Response Times¶