import numpy as np

# Example: Two dice rolls
# P(sum = 7) via simulation (law of large numbers)
np.random.seed(42)
n_trials = 100_000
die1 = np.random.randint(1, 7, n_trials)
die2 = np.random.randint(1, 7, n_trials)
sums = die1 + die2

p_sum_7 = (sums == 7).mean()
print(f"P(sum=7) simulated:     {p_sum_7:.4f}")
print(f"P(sum=7) theoretical:   {6/36:.4f}")   # 6 ways out of 36 total

# P(at least one 6) = 1 - P(no 6 on either die)
p_no_six = (die1 != 6).mean() * (die2 != 6).mean()
p_at_least_one_six = 1 - p_no_six
print(f"P(at least one 6):      {p_at_least_one_six:.4f}")
print(f"Theoretical (1 - 25/36):{1 - 25/36:.4f}")

# Union: P(sum=7 OR sum=11)
p_7_or_11 = ((sums == 7) | (sums == 11)).mean()
print(f"P(sum=7 or sum=11):     {p_7_or_11:.4f}")
print(f"Theoretical (8/36):     {8/36:.4f}")

# Bayes' theorem for medical testing
prevalence   = 0.01   # P(Disease) = 1% of population
sensitivity  = 0.95   # P(Positive | Disease)
specificity  = 0.90   # P(Negative | No Disease) → P(Positive | No Disease) = 0.10

# ── Manual Bayes calculation ──
# P(B) = P(Positive) = P(Pos|Disease)*P(Disease) + P(Pos|No Disease)*P(No Disease)
p_positive = sensitivity * prevalence + (1 - specificity) * (1 - prevalence)

# P(Disease | Positive) = P(Positive | Disease) * P(Disease) / P(Positive)
p_disease_given_positive = (sensitivity * prevalence) / p_positive

print(f"P(Positive):              {p_positive:.4f}  ({p_positive:.1%})")
print(f"P(Disease | Test+):       {p_disease_given_positive:.4f}  ({p_disease_given_positive:.1%})")
# Result: ~8.7% — even with a positive test, most positives are FALSE POSITIVES

# ── Intuition via 10,000 people ──
n = 10_000
diseased     = int(n * prevalence)           # 100 people with disease
healthy      = n - diseased                  # 9900 people healthy

true_positive  = int(diseased * sensitivity)           # 95 correctly detected
false_positive = int(healthy * (1 - specificity))      # 990 healthy but test positive
all_positives  = true_positive + false_positive        # 1085 total positives

ppv = true_positive / all_positives   # Positive Predictive Value
print(f"\nOut of 10,000 people tested:")
print(f"  True positives:  {true_positive}")
print(f"  False positives: {false_positive}")
print(f"  P(Disease | Positive): {ppv:.1%}")   # ~8.8%
print("\nThis is the base rate fallacy: ignoring the rarity of disease inflates perceived accuracy")

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

# ── Binomial: number of defects in a batch of 100 items ──
n = 100      # batch size
p = 0.03     # defect probability per item (3%)

binom = stats.binom(n=n, p=p)

# PMF: probability of exactly k defects
k_values = np.arange(0, 16)
pmf = binom.pmf(k_values)

print(f"P(0 defects):   {binom.pmf(0):.4f}")
print(f"P(exactly 3):   {binom.pmf(3):.4f}")
print(f"P(at most 5):   {binom.cdf(5):.4f}")   # CDF: P(X ≤ 5)
print(f"Expected defects (mean): {binom.mean():.1f}")  # n*p = 3

# ── Plot PMF ──
fig, ax = plt.subplots(figsize=(9, 4))
ax.bar(k_values, pmf, color='steelblue', edgecolor='white', alpha=0.8)
ax.set_title('Binomial(n=100, p=0.03) — Defects per Batch', fontsize=12)
ax.set_xlabel('Number of Defects'); ax.set_ylabel('Probability')
ax.set_xlim(-0.5, 15.5)
plt.tight_layout(); plt.show()

from scipy import stats
import numpy as np

# ── Poisson: number of server errors per hour, historical average λ=4 ──
lam = 4
poisson = stats.poisson(mu=lam)

print(f"P(exactly 0 errors):     {poisson.pmf(0):.4f}")
print(f"P(exactly 4 errors):     {poisson.pmf(4):.4f}")   # highest single-value probability
print(f"P(more than 8 errors):   {1 - poisson.cdf(8):.4f}")   # tail risk — worth alerting on?
print(f"Mean = Variance = λ:     mean={poisson.mean():.1f}, var={poisson.var():.1f}")

# Poisson as the n->infinity, p->0 limit of Binomial with n*p held constant
n_large, p_small = 10_000, lam / 10_000
binom_approx = stats.binom(n=n_large, p=p_small)
print(f"\nBinomial({n_large}, {p_small:.5f}) P(X=4): {binom_approx.pmf(4):.4f}  (should ≈ Poisson's {poisson.pmf(4):.4f})")

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

# ── Uniform distribution ──
uniform_rv = stats.uniform(loc=0, scale=10)   # Uniform[0, 10]
print(f"P(2 ≤ X ≤ 5) = {uniform_rv.cdf(5) - uniform_rv.cdf(2):.2f}")  # 0.30
print(f"E[X] = {uniform_rv.mean():.1f}")   # 5.0

# ── Normal distribution ──
normal_rv = stats.norm(loc=70, scale=10)   # mean=70, std=10
print(f"P(X ≤ 85) = {normal_rv.cdf(85):.4f}")    # ~0.9332
print(f"P(60 ≤ X ≤ 80) = {normal_rv.cdf(80) - normal_rv.cdf(60):.4f}")  # ~0.6827

# ── Plot both ──
x_unif = np.linspace(-1, 11, 300)
x_norm = np.linspace(30, 110, 300)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 4))
ax1.plot(x_unif, uniform_rv.pdf(x_unif), 'steelblue', lw=2)
ax1.set_title('Uniform(0, 10) PDF'); ax1.set_ylim(0, 0.15)

ax2.plot(x_norm, normal_rv.pdf(x_norm), 'steelblue', lw=2)
ax2.set_title('Normal(μ=70, σ=10) PDF')
plt.tight_layout(); plt.show()

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

np.random.seed(0)
symmetric    = np.random.normal(50, 10, 2000)
right_skewed = np.random.exponential(scale=20, size=2000)   # income-like
left_skewed  = 100 - np.random.exponential(scale=8, size=2000)

datasets = {
    'Symmetric (Normal)':     symmetric,
    'Right-Skewed (Positive)': right_skewed,
    'Left-Skewed (Negative)':  left_skewed
}

for name, data in datasets.items():
    sk = pd.Series(data).skew()
    mean = np.mean(data)
    median = np.median(data)
    print(f"{name:30s}  skew={sk:+.2f}  mean={mean:.1f}  median={median:.1f}")

# Symmetric (Normal)              skew=-0.04  mean=50.0  median=50.1
# Right-Skewed (Positive)         skew=+2.05  mean=20.5  median=13.9  (mean > median)
# Left-Skewed (Negative)          skew=-2.30  mean=87.9  median=94.3  (mean < median)

# ── Rule of thumb: when to apply log transform ──
print(f"\nRight-skewed skew: {pd.Series(right_skewed).skew():.2f}")
log_transformed = np.log1p(right_skewed)
print(f"After log1p:       {pd.Series(log_transformed).skew():.2f}")   # near 0

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

# ── Normal distribution with scipy ──
mu, sigma = 170, 8   # human height: mean 170cm, std 8cm
normal = stats.norm(loc=mu, scale=sigma)

# PDF: probability density at specific points
x = np.linspace(140, 200, 300)
pdf_values = normal.pdf(x)

# CDF: probability that X ≤ value
p_below_180 = normal.cdf(180)
print(f"P(height ≤ 180cm) = {p_below_180:.4f} ({p_below_180:.1%})")

# P(a ≤ X ≤ b) = CDF(b) - CDF(a)
p_between  = normal.cdf(178) - normal.cdf(162)
print(f"P(162 ≤ height ≤ 178cm) = {p_between:.4f} ({p_between:.1%})")

# Inverse CDF (PPF — percent point function): what value has exactly p% below it?
p90_height = normal.ppf(0.90)
print(f"90th percentile height = {p90_height:.1f} cm")

# ── 68-95-99.7 rule verification ──
for n_std in [1, 2, 3]:
    lo = normal.cdf(mu - n_std * sigma)
    hi = normal.cdf(mu + n_std * sigma)
    p = hi - lo
    print(f"μ ± {n_std}σ = [{mu - n_std*sigma}, {mu + n_std*sigma}] cm — {p:.4f} ({p:.1%})")

# ── Plot ──
fig, ax = plt.subplots(figsize=(10, 5))
ax.plot(x, pdf_values, 'steelblue', lw=2.5, label='Normal(170, 8)')
ax.fill_between(x, pdf_values, where=(x >= mu - sigma) & (x <= mu + sigma),
                alpha=0.3, color='steelblue', label='68% (±1σ)')
ax.fill_between(x, pdf_values, where=((x >= mu - 2*sigma) & (x < mu - sigma)) |
                                     ((x > mu + sigma) & (x <= mu + 2*sigma)),
                alpha=0.2, color='orange', label='95% (±2σ)')
ax.axvline(mu, color='red', linestyle='--', linewidth=1.5, label=f'μ = {mu}')
ax.set_title('Normal Distribution — 68-95-99.7 Rule', fontsize=12, fontweight='bold')
ax.set_xlabel('Height (cm)'); ax.set_ylabel('Probability Density')
ax.legend(); plt.tight_layout(); plt.show()

from scipy import stats
import numpy as np

# Standard normal
std_normal = stats.norm(loc=0, scale=1)

# Look up probabilities from the standard normal table
print(f"P(Z ≤ 0)    = {std_normal.cdf(0):.4f}")   # 0.5000 (symmetric at 0)
print(f"P(Z ≤ 1.96) = {std_normal.cdf(1.96):.4f}")  # 0.9750 (95% CI upper bound)
print(f"P(Z ≤ 2)    = {std_normal.cdf(2):.4f}")   # 0.9772

# Converting from a regular normal to standard normal
mu, sigma = 170, 8   # heights
x_val = 182          # a specific height

z_score = (x_val - mu) / sigma
print(f"\nx = {x_val} cm → z = ({x_val} − {mu}) / {sigma} = {z_score:.2f}")
print(f"P(height ≤ {x_val}) = {stats.norm.cdf(z_score):.4f}")

# Both approaches give the same answer:
print(f"Directly:   P(height ≤ {x_val}) = {stats.norm.cdf(x_val, loc=mu, scale=sigma):.4f}")

# Critical z-values you'll encounter constantly:
print("\nCritical z-values:")
for p in [0.90, 0.95, 0.975, 0.99]:
    z = stats.norm.ppf(p)
    print(f"  z for P(Z ≤ z) = {p:.3f}:  z = {z:.3f}")

import numpy as np
import pandas as pd

np.random.seed(42)
# Simulate a dataset with a few outliers
data = np.append(np.random.normal(loc=100, scale=15, size=200),
                 [200, 210, -50])   # outliers

df = pd.DataFrame({'value': data})

# ── Compute Z-scores manually ──
mean  = df['value'].mean()
std   = df['value'].std()
df['z_score'] = (df['value'] - mean) / std

print(f"Mean:  {mean:.2f}")
print(f"Std:   {std:.2f}")
print(f"\nTop 5 highest z-scores:")
print(df.nlargest(5, 'z_score')[['value', 'z_score']].round(2))

# ── Flag outliers using |z| > 3 threshold ──
df['is_outlier'] = df['z_score'].abs() > 3
print(f"\nOutliers detected: {df['is_outlier'].sum()}")
print(df[df['is_outlier']][['value', 'z_score']].round(2))

# ── Using scipy for standardization ──
from scipy import stats
z_scores_scipy = stats.zscore(data)
print(f"\nFirst 5 z-scores (scipy): {z_scores_scipy[:5].round(3)}")

# ── Or do it yourself with NumPy (same math, no library) ──
X = df['value'].values
X_scaled = (X - X.mean()) / X.std()

print(f"\nAfter standardization:")
print(f"  Mean:  {X_scaled.mean():.6f}")   # ~0.0
print(f"  Std:   {X_scaled.std():.6f}")    # ~1.0
# In Phase 2 (Lesson 10) you'll meet scikit-learn's StandardScaler,
# which wraps exactly this computation in a reusable tool

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

np.random.seed(0)
# Right-skewed data: exponential with outliers
data = np.append(np.random.exponential(scale=50, size=300),
                 [600, 700, 750])   # known outliers

df = pd.DataFrame({'value': data})

# ── Method 1: Z-score ──
mean, std = df['value'].mean(), df['value'].std()
df['z_score']      = (df['value'] - mean) / std
df['z_outlier']    = df['z_score'].abs() > 3

# ── Method 2: IQR ──
q1, q3 = df['value'].quantile(0.25), df['value'].quantile(0.75)
iqr = q3 - q1
df['iqr_outlier'] = (df['value'] < q1 - 1.5*iqr) | (df['value'] > q3 + 1.5*iqr)

print(f"Z-score  outliers detected: {df['z_outlier'].sum()}")
print(f"IQR      outliers detected: {df['iqr_outlier'].sum()}")

# On skewed data, IQR tends to flag more true outliers from the right tail
# Z-score may miss them if the extreme values inflate the std

# ── Compare which known outliers each method catches ──
for method, col in [('Z-score', 'z_outlier'), ('IQR', 'iqr_outlier')]:
    caught = df[df['value'] > 500][col].sum()
    print(f"{method}: caught {caught}/3 extreme outliers (>500)")

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)

def pd_skew(x):
    import pandas as pd
    return pd.Series(x).skew()

# A deliberately NON-normal population: exponential (heavily right-skewed)
population = np.random.exponential(scale=2.0, size=100_000)
print(f"Population mean={population.mean():.3f}, skewness={pd_skew(population):.2f} "
      f"(heavily right-skewed, nothing like a bell curve)")

# Draw many samples of size n, record each sample's MEAN
def sampling_distribution_of_the_mean(population, n, n_samples=5000):
    means = [np.random.choice(population, size=n, replace=True).mean()
             for _ in range(n_samples)]
    return np.array(means)

fig, axes = plt.subplots(1, 4, figsize=(16, 3.5))
axes[0].hist(population, bins=60, color='indianred', alpha=0.8)
axes[0].set_title(f'Population\n(skew={pd_skew(population):.2f}, clearly NOT normal)')

for ax, n in zip(axes[1:], [2, 10, 50]):
    sample_means = sampling_distribution_of_the_mean(population, n)
    ax.hist(sample_means, bins=60, color='steelblue', alpha=0.8)
    ax.set_title(f'Means of {n}-sample draws\n(skew={pd_skew(sample_means):.2f})')

plt.tight_layout()
plt.savefig('clt_demo.png', dpi=150)
# As n grows from 2 -> 10 -> 50, the histogram of sample means becomes
# visibly more symmetric and bell-shaped, even though every single sample
# was drawn from a badly skewed exponential population.