Week 6: Probability, Distributions and Statistical Inference

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

# Basic probability
# Probability = Number of favorable outcomes / Total possible outcomes

print("="*50)
print("BASIC PROBABILITY EXAMPLES")
print("="*50)

# Example 1: Rolling a die
# Probability of rolling a 4
favorable_outcomes = 1
total_outcomes = 6
prob_four = favorable_outcomes / total_outcomes
print(f"\nP(rolling a 4) = {favorable_outcomes}/{total_outcomes} = {prob_four:.3f} or {prob_four*100:.1f}%")

# Example 2: Drawing from a deck
# Probability of drawing a red card
red_cards = 26
total_cards = 52
prob_red = red_cards / total_cards
print(f"P(red card) = {red_cards}/{total_cards} = {prob_red:.3f} or {prob_red*100:.1f}%")

# Example 3: Probability rules
print("\n" + "="*50)
print("PROBABILITY RULES")
print("="*50)

# Addition rule: P(A or B) = P(A) + P(B) - P(A and B)
P_A = 0.3 # Probability of event A
P_B = 0.4 # Probability of event B
P_A_and_B = 0.1 # Probability of both occurring

P_A_or_B = P_A + P_B - P_A_and_B
print(f"\nAddition Rule: P(A or B) = P(A) + P(B) - P(A and B)")
print(f"P(A or B) = {P_A} + {P_B} - {P_A_and_B} = {P_A_or_B}")

# Multiplication rule: P(A and B) = P(A) * P(B|A)
# For independent events: P(A and B) = P(A) * P(B)
P_heads = 0.5
P_tails = 0.5
P_two_heads = P_heads * P_heads
print(f"\nMultiplication Rule (independent): P(two heads) = 0.5 × 0.5 = {P_two_heads}")

# Complement rule: P(not A) = 1 - P(A)
P_rain = 0.3
P_no_rain = 1 - P_rain
print(f"\nComplement Rule: P(no rain) = 1 - P(rain) = 1 - {P_rain} = {P_no_rain}")

# Conditional probability: P(A|B) = P(A and B) / P(B)
print("\n" + "="*50)
print("CONDITIONAL PROBABILITY")
print("="*50)

# Example: Disease test accuracy
# P(Disease|Positive Test)
P_positive_given_disease = 0.95 # Sensitivity
P_negative_given_no_disease = 0.90 # Specificity
P_disease = 0.01 # Prior probability of disease

# Using Bayes' theorem
P_positive = (P_positive_given_disease * P_disease) + ((1 - P_negative_given_no_disease) * (1 - P_disease))
P_disease_given_positive = (P_positive_given_disease * P_disease) / P_positive

print(f"Sensitivity (P(+|disease)): {P_positive_given_disease}")
print(f"Specificity (P(-|no disease)): {P_negative_given_no_disease}")
print(f"P(disease) prior: {P_disease}")
print(f"P(+) overall: {P_positive:.4f}")
print(f"P(disease|+) posterior: {P_disease_given_positive:.4f}")
print(f"\nInterpretation: Even with positive test, only {P_disease_given_positive*100:.1f}% probability of having disease!")

# Simulation: Rolling dice
print("\n" + "="*50)
print("SIMULATION: ROLLING DICE")
print("="*50)

np.random.seed(42)
num_rolls = 100000

# Simulate rolling a die
rolls = np.random.randint(1, 7, num_rolls)

# Calculate probabilities
for face in range(1, 7):
prob = (rolls == face).sum() / num_rolls
print(f"P(rolling {face}) = {prob:.4f}")

print(f"\nExpected: {1/6:.4f} for each face")

# Simulate coin flips
print("\n" + "="*50)
print("SIMULATION: COIN FLIPS")
print("="*50)

flips = np.random.randint(0, 2, 100000)
prob_heads = (flips == 1).sum() / 100000
prob_tails = (flips == 0).sum() / 100000

print(f"P(Heads) = {prob_heads:.4f}")
print(f"P(Tails) = {prob_tails:.4f}")

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

print("="*50)
print("DISCRETE DISTRIBUTIONS")
print("="*50)

# 1. BERNOULLI DISTRIBUTION
# Single trial with two outcomes (success/failure)
print("\n1. BERNOULLI DISTRIBUTION")
print("Example: Single coin flip, single product working/failing")

p = 0.7 # Probability of success
bernoulli = stats.bernoulli(p)

print(f"P(success) = {p}")
print(f"P(failure) = {1-p}")
print(f"Mean = {bernoulli.mean()}")
print(f"Variance = {bernoulli.var()}")

# 2. BINOMIAL DISTRIBUTION
# Number of successes in n independent trials
print("\n" + "="*50)
print("2. BINOMIAL DISTRIBUTION")
print("Example: Number of heads in 10 coin flips, defects in 100 items")

n = 10 # Number of trials
p = 0.5 # Probability of success per trial
binomial = stats.binom(n, p)

print(f"\nParameters: n={n} trials, p={p}")
print(f"Mean = {binomial.mean()}")
print(f"Variance = {binomial.var()}")

# Probability of exactly 5 successes
prob_5 = binomial.pmf(5)
print(f"P(X=5) = {prob_5:.4f}")

# Probability of 5 or fewer successes
prob_le5 = binomial.cdf(5)
print(f"P(X≤5) = {prob_le5:.4f}")

# Simulate
simulated = np.random.binomial(n, p, 10000)
print(f"\nSimulated mean: {simulated.mean():.3f}")
print(f"Actual mean: {binomial.mean():.3f}")

# 3. POISSON DISTRIBUTION
# Number of events in fixed interval
print("\n" + "="*50)
print("3. POISSON DISTRIBUTION")
print("Example: Number of customer arrivals per hour, defects per batch")

lambda_param = 3 # Average number of events
poisson = stats.poisson(lambda_param)

print(f"\nParameter (λ) = {lambda_param}")
print(f"Mean = {poisson.mean()}")
print(f"Variance = {poisson.var()}")

# Probability of exactly 5 events
prob_5 = poisson.pmf(5)
print(f"P(X=5) = {prob_5:.4f}")

# Probability of 5 or fewer events
prob_le5 = poisson.cdf(5)
print(f"P(X≤5) = {prob_le5:.4f}")

# Simulate
simulated = np.random.poisson(lambda_param, 10000)
print(f"\nSimulated mean: {simulated.mean():.3f}")
print(f"Actual mean: {poisson.mean():.3f}")

# Visualization
fig, axes = plt.subplots(1, 3, figsize=(15, 4))

# Binomial
x = range(0, n+1)
pmf = binomial.pmf(x)
axes[0].bar(x, pmf, alpha=0.7, color='blue')
axes[0].set_title(f'Binomial(n={n}, p={p})')
axes[0].set_xlabel('Number of successes')
axes[0].set_ylabel('Probability')
axes[0].grid(True, alpha=0.3)

# Poisson
x = range(0, 10)
pmf = poisson.pmf(x)
axes[1].bar(x, pmf, alpha=0.7, color='green')
axes[1].set_title(f'Poisson(λ={lambda_param})')
axes[1].set_xlabel('Number of events')
axes[1].set_ylabel('Probability')
axes[1].grid(True, alpha=0.3)

# Distribution comparison
axes[2].bar([0, 1], [1-p, p], alpha=0.7, color='red', label='Bernoulli')
axes[2].set_title(f'Bernoulli(p={p})')
axes[2].set_xlabel('Outcome')
axes[2].set_ylabel('Probability')
axes[2].set_xticks([0, 1])
axes[2].set_xticklabels(['Failure', 'Success'])
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Real-world example: Product quality
print("\n" + "="*50)
print("REAL-WORLD: PRODUCT QUALITY CONTROL")
print("="*50)

# Defect rate is 2%, inspecting 50 units
n_units = 50
defect_rate = 0.02
quality = stats.binom(n_units, defect_rate)

print(f"\nDefect rate: {defect_rate*100}%")
print(f"Inspecting: {n_units} units")
print(f"Expected defects: {quality.mean()}")
print(f"P(0 defects) = {quality.pmf(0):.4f}")
print(f"P(≥2 defects) = {1 - quality.cdf(1):.4f}")

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

print("="*50)
print("CONTINUOUS DISTRIBUTIONS")
print("="*50)

# 1. NORMAL (GAUSSIAN) DISTRIBUTION
print("\n1. NORMAL DISTRIBUTION")
print("Example: Heights, test scores, measurement errors")

mu = 100 # Mean
sigma = 15 # Standard deviation
normal = stats.norm(mu, sigma)

print(f"\nParameters: μ={mu}, σ={sigma}")
print(f"Mean = {normal.mean()}")
print(f"Variance = {normal.var()}")

# Probabilities
prob_less_100 = normal.cdf(100)
print(f"P(X < 100) = {prob_less_100:.4f}")

prob_100_110 = normal.cdf(110) - normal.cdf(100)
print(f"P(100 < X < 110) = {prob_100_110:.4f}")

# Percentiles
percentile_90 = normal.ppf(0.90)
print(f"90th percentile = {percentile_90:.2f}")

# Z-scores (standardized)
z_score = (110 - mu) / sigma
print(f"Z-score for X=110: {z_score:.3f}")

# Standard normal (μ=0, σ=1)
print(f"\nP(Z < 0) in standard normal = {stats.norm.cdf(0):.4f}")
print(f"P(-1 < Z < 1) = {stats.norm.cdf(1) - stats.norm.cdf(-1):.4f}")
print(f"P(-2 < Z < 2) = {stats.norm.cdf(2) - stats.norm.cdf(-2):.4f}")
print(f"P(-3 < Z < 3) = {stats.norm.cdf(3) - stats.norm.cdf(-3):.4f}")

# 2. EXPONENTIAL DISTRIBUTION
print("\n" + "="*50)
print("2. EXPONENTIAL DISTRIBUTION")
print("Example: Time until next event, equipment failure time")

lambda_param = 0.5 # Rate parameter
exponential = stats.expon(scale=1/lambda_param)

print(f"\nParameter (λ) = {lambda_param}")
print(f"Mean = {exponential.mean():.3f}")
print(f"Variance = {exponential.var():.3f}")

# Probability density
x = np.linspace(0, 10, 100)
pdf = exponential.pdf(x)

# Cumulative probability
prob_less_2 = exponential.cdf(2)
print(f"P(X < 2) = {prob_less_2:.4f}")

# 3. UNIFORM DISTRIBUTION
print("\n" + "="*50)
print("3. UNIFORM DISTRIBUTION")
print("Example: Random number generation, continuous random selection")

a, b = 0, 10 # Bounds
uniform = stats.uniform(a, b-a)

print(f"\nBounds: [{a}, {b}]")
print(f"Mean = {uniform.mean():.3f}")
print(f"Variance = {uniform.var():.3f}")

prob_3_7 = uniform.cdf(7) - uniform.cdf(3)
print(f"P(3 < X < 7) = {prob_3_7:.4f}")

# Visualization
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# Normal distribution
x = np.linspace(mu - 4*sigma, mu + 4*sigma, 100)
pdf = normal.pdf(x)
axes[0, 0].plot(x, pdf, 'b-', linewidth=2, label='Normal')
axes[0, 0].fill_between(x, pdf, alpha=0.3)
axes[0, 0].axvline(mu, color='r', linestyle='--', label=f'Mean = {mu}')
axes[0, 0].set_title(f'Normal Distribution μ={mu}, σ={sigma}')
axes[0, 0].set_xlabel('Value')
axes[0, 0].set_ylabel('Probability Density')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

# Exponential distribution
x = np.linspace(0, 10, 100)
pdf = exponential.pdf(x)
axes[0, 1].plot(x, pdf, 'g-', linewidth=2)
axes[0, 1].fill_between(x, pdf, alpha=0.3, color='green')
axes[0, 1].set_title(f'Exponential Distribution λ={lambda_param}')
axes[0, 1].set_xlabel('Value')
axes[0, 1].set_ylabel('Probability Density')
axes[0, 1].grid(True, alpha=0.3)

# Uniform distribution
x = np.linspace(a-1, b+1, 100)
pdf = uniform.pdf(x)
axes[1, 0].plot(x, pdf, 'r-', linewidth=2)
axes[1, 0].fill_between(x, pdf, alpha=0.3, color='red')
axes[1, 0].set_title(f'Uniform Distribution [{a}, {b}]')
axes[1, 0].set_xlabel('Value')
axes[1, 0].set_ylabel('Probability Density')
axes[1, 0].grid(True, alpha=0.3)

# Multiple normal distributions for comparison
x = np.linspace(-50, 250, 100)
for mu_i in [100, 150, 200]:
pdf = stats.norm.pdf(x, mu_i, 15)
axes[1, 1].plot(x, pdf, linewidth=2, label=f'μ={mu_i}')
axes[1, 1].set_title('Normal Distributions with Different Means')
axes[1, 1].set_xlabel('Value')
axes[1, 1].set_ylabel('Probability Density')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Real-world example: Test scores
print("\n" + "="*50)
print("REAL-WORLD: TEST SCORE ANALYSIS")
print("="*50)

mean_score = 75
std_score = 8
test_dist = stats.norm(mean_score, std_score)

print(f"\nMean score: {mean_score}")
print(f"Standard deviation: {std_score}")
print(f"P(score < 70) = {test_dist.cdf(70):.4f} or {test_dist.cdf(70)*100:.2f}%")
print(f"P(score > 90) = {1 - test_dist.cdf(90):.4f} or {(1-test_dist.cdf(90))*100:.2f}%")
print(f"P(70 < score < 85) = {test_dist.cdf(85) - test_dist.cdf(70):.4f} or {(test_dist.cdf(85) - test_dist.cdf(70))*100:.2f}%")

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

print("="*50)
print("SAMPLING DISTRIBUTIONS")
print("="*50)

print("""
KEY CONCEPTS:
- Population: All items of interest
- Sample: Subset of population
- Sampling distribution: Distribution of sample statistics
- Standard Error: Standard deviation of sampling distribution
""")

# Create a population
np.random.seed(42)
population = np.random.normal(100, 15, 100000)

print(f"\nPopulation statistics:")
print(f"Mean: {population.mean():.2f}")
print(f"Std Dev: {population.std():.2f}")

# Draw many samples and calculate means
sample_size = 30
num_samples = 10000
sample_means = []

for i in range(num_samples):
sample = np.random.choice(population, sample_size, replace=False)
sample_means.append(sample.mean())

sample_means = np.array(sample_means)

print(f"\n" + "="*50)
print("SAMPLING DISTRIBUTION OF MEANS")
print("="*50)
print(f"Sample size: {sample_size}")
print(f"Number of samples: {num_samples}")
print(f"\nMean of sample means: {sample_means.mean():.2f}")
print(f"Std Dev of sample means (Standard Error): {sample_means.std():.2f}")

# Compare to theoretical standard error
se_theoretical = population.std() / np.sqrt(sample_size)
print(f"Theoretical SE: {se_theoretical:.2f}")

# Effect of sample size
print("\n" + "="*50)
print("EFFECT OF SAMPLE SIZE")
print("="*50)

sample_sizes = [10, 30, 100, 300]
results = {}

for n in sample_sizes:
sample_means_n = []
for i in range(5000):
sample = np.random.choice(population, n, replace=False)
sample_means_n.append(sample.mean())

results[n] = np.array(sample_means_n)
se = population.std() / np.sqrt(n)
print(f"Sample size {n:3d}: SE = {results[n].std():.2f} (theoretical: {se:.2f})")

# Visualization
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

for idx, n in enumerate(sample_sizes):
row = idx // 2
col = idx % 2

axes[row, col].hist(results[n], bins=50, alpha=0.7, color='blue', edgecolor='black')
axes[row, col].set_title(f'Sample Size = {n}')
axes[row, col].set_xlabel('Sample Mean')
axes[row, col].set_ylabel('Frequency')
axes[row, col].axvline(population.mean(), color='red', linestyle='--', label='Population Mean')
axes[row, col].legend()
axes[row, col].grid(True, alpha=0.3)

plt.suptitle('Sampling Distribution of Means (Different Sample Sizes)', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()

# Standard Error for different statistics
print("\n" + "="*50)
print("STANDARD ERROR FORMULAS")
print("="*50)

print("""
For sample mean: SE = σ / √n
For sample proportion: SE = √(p(1-p) / n)
For sample variance: More complex formula
""")

# Example: Proportion
p = 0.4 # True population proportion
n = 100
se_proportion = np.sqrt(p * (1-p) / n)
print(f"\nFor proportion: p={p}, n={n}")
print(f"SE = √({p}×{1-p}/{n}) = {se_proportion:.4f}")

# Simulate
sample_proportions = []
for i in range(10000):
sample = np.random.binomial(n, p) / n
sample_proportions.append(sample)

print(f"Simulated SE: {np.std(sample_proportions):.4f}")
print(f"Theoretical SE: {se_proportion:.4f}")

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

print("="*50)
print("CENTRAL LIMIT THEOREM (CLT)")
print("="*50)

print("""
THE CENTRAL LIMIT THEOREM STATES:
If you take repeated samples from ANY distribution,
the DISTRIBUTION OF THE SAMPLE MEANS approaches
a NORMAL DISTRIBUTION as sample size increases.

This is true regardless of the underlying distribution!

Implications:
- We can make inferences about populations using sample means
- Normal distribution can approximate many distributions
- Sample size matters: larger samples → more normal
""")

# Demonstrate CLT with different distributions
fig, axes = plt.subplots(3, 3, figsize=(15, 12))

# 1. UNIFORM DISTRIBUTION
print("\n" + "="*50)
print("1. UNIFORM DISTRIBUTION")
print("="*50)

# Generate population (uniform)
population_uniform = np.random.uniform(0, 10, 50000)
axes[0, 0].hist(population_uniform, bins=50, alpha=0.7, color='blue', edgecolor='black')
axes[0, 0].set_title('Population (Uniform)')
axes[0, 0].set_xlabel('Value')
axes[0, 0].set_ylabel('Frequency')
axes[0, 0].grid(True, alpha=0.3)

# Sample means with n=10
sample_means_10 = []
for i in range(5000):
sample = np.random.choice(population_uniform, 10)
sample_means_10.append(sample.mean())

axes[0, 1].hist(sample_means_10, bins=50, alpha=0.7, color='green', edgecolor='black')
axes[0, 1].set_title('Sample Means (n=10)')
axes[0, 1].set_xlabel('Sample Mean')
axes[0, 1].set_ylabel('Frequency')
axes[0, 1].grid(True, alpha=0.3)

# Sample means with n=100
sample_means_100 = []
for i in range(5000):
sample = np.random.choice(population_uniform, 100)
sample_means_100.append(sample.mean())

axes[0, 2].hist(sample_means_100, bins=50, alpha=0.7, color='red', edgecolor='black')
axes[0, 2].set_title('Sample Means (n=100)')
axes[0, 2].set_xlabel('Sample Mean')
axes[0, 2].set_ylabel('Frequency')
axes[0, 2].grid(True, alpha=0.3)

# 2. EXPONENTIAL DISTRIBUTION (Right-skewed)
print("\n2. EXPONENTIAL DISTRIBUTION")

population_exp = np.random.exponential(2, 50000)
axes[1, 0].hist(population_exp, bins=50, alpha=0.7, color='blue', edgecolor='black')
axes[1, 0].set_title('Population (Exponential)')
axes[1, 0].set_xlabel('Value')
axes[1, 0].set_ylabel('Frequency')
axes[1, 0].grid(True, alpha=0.3)

sample_means_exp_10 = []
for i in range(5000):
sample = np.random.choice(population_exp, 10)
sample_means_exp_10.append(sample.mean())

axes[1, 1].hist(sample_means_exp_10, bins=50, alpha=0.7, color='green', edgecolor='black')
axes[1, 1].set_title('Sample Means (n=10)')
axes[1, 1].set_xlabel('Sample Mean')
axes[1, 1].set_ylabel('Frequency')
axes[1, 1].grid(True, alpha=0.3)

sample_means_exp_100 = []
for i in range(5000):
sample = np.random.choice(population_exp, 100)
sample_means_exp_100.append(sample.mean())

axes[1, 2].hist(sample_means_exp_100, bins=50, alpha=0.7, color='red', edgecolor='black')
axes[1, 2].set_title('Sample Means (n=100)')
axes[1, 2].set_xlabel('Sample Mean')
axes[1, 2].set_ylabel('Frequency')
axes[1, 2].grid(True, alpha=0.3)

# 3. BIMODAL DISTRIBUTION
print("\n3. BIMODAL DISTRIBUTION")

population_bimodal = np.concatenate([
np.random.normal(20, 2, 25000),
np.random.normal(80, 2, 25000)
])

axes[2, 0].hist(population_bimodal, bins=50, alpha=0.7, color='blue', edgecolor='black')
axes[2, 0].set_title('Population (Bimodal)')
axes[2, 0].set_xlabel('Value')
axes[2, 0].set_ylabel('Frequency')
axes[2, 0].grid(True, alpha=0.3)

sample_means_bi_10 = []
for i in range(5000):
sample = np.random.choice(population_bimodal, 10)
sample_means_bi_10.append(sample.mean())

axes[2, 1].hist(sample_means_bi_10, bins=50, alpha=0.7, color='green', edgecolor='black')
axes[2, 1].set_title('Sample Means (n=10)')
axes[2, 1].set_xlabel('Sample Mean')
axes[2, 1].set_ylabel('Frequency')
axes[2, 1].grid(True, alpha=0.3)

sample_means_bi_100 = []
for i in range(5000):
sample = np.random.choice(population_bimodal, 100)
sample_means_bi_100.append(sample.mean())

axes[2, 2].hist(sample_means_bi_100, bins=50, alpha=0.7, color='red', edgecolor='black')
axes[2, 2].set_title('Sample Means (n=100)')
axes[2, 2].set_xlabel('Sample Mean')
axes[2, 2].set_ylabel('Frequency')
axes[2, 2].grid(True, alpha=0.3)

plt.suptitle('Central Limit Theorem: Sample Means Approach Normal Distribution',
fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()

print("\n" + "="*50)
print("OBSERVATIONS")
print("="*50)
print("As sample size increases, sample means distribution:")
print("1. Becomes more bell-shaped (normal)")
print("2. Becomes more concentrated (smaller variance)")
print("3. Approaches normal regardless of population shape")

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

print("="*50)
print("CONFIDENCE INTERVALS")
print("="*50)

print("""
WHAT IS A CONFIDENCE INTERVAL?
- A range of values estimated to contain the true population parameter
- Constructed from sample data
- Associated with a confidence level (e.g., 95%)

INTERPRETATION:
If we repeated our sampling 100 times and created CI each time,
95% of those intervals would contain the true population parameter.

Common confidence levels: 90%, 95%, 99%
""")

# Example 1: Confidence interval for mean
print("\nEXAMPLE 1: CI FOR POPULATION MEAN")
print("="*50)

# Sample data
sample_data = np.array([85, 92, 78, 88, 95, 82, 90, 87, 93, 80])
print(f"Sample: {sample_data}")
print(f"n = {len(sample_data)}")

# Calculate statistics
sample_mean = sample_data.mean()
sample_std = sample_data.std(ddof=1) # Sample std dev
se = sample_std / np.sqrt(len(sample_data))

print(f"Sample mean: {sample_mean:.2f}")
print(f"Sample std dev: {sample_std:.2f}")
print(f"Standard error: {se:.2f}")

# 95% Confidence interval using t-distribution
confidence_level = 0.95
alpha = 1 - confidence_level
t_critical = stats.t.ppf(1 - alpha/2, df=len(sample_data)-1)
margin_error = t_critical * se

ci_lower = sample_mean - margin_error
ci_upper = sample_mean + margin_error

print(f"\nt-critical (df={len(sample_data)-1}): {t_critical:.3f}")
print(f"Margin of error: {margin_error:.2f}")
print(f"95% CI: [{ci_lower:.2f}, {ci_upper:.2f}]")

# Using scipy
ci = stats.t.interval(0.95, len(sample_data)-1, loc=sample_mean, scale=se)
print(f"Using scipy: [{ci[0]:.2f}, {ci[1]:.2f}]")

# Example 2: Confidence interval for proportion
print("\n" + "="*50)
print("EXAMPLE 2: CI FOR PROPORTION")
print("="*50)

# Survey: 200 people, 120 prefer product A
n = 200
x = 120 # Number preferring A
p_hat = x / n

print(f"Sample size: {n}")
print(f"Number preferring: {x}")
print(f"Sample proportion: {p_hat:.3f}")

# 95% CI for proportion
se_prop = np.sqrt(p_hat * (1 - p_hat) / n)
z_critical = stats.norm.ppf(0.975) # For 95%
margin_error_prop = z_critical * se_prop

ci_lower_prop = p_hat - margin_error_prop
ci_upper_prop = p_hat + margin_error_prop

print(f"Standard error: {se_prop:.4f}")
print(f"z-critical: {z_critical:.3f}")
print(f"Margin of error: {margin_error_prop:.3f}")
print(f"95% CI: [{ci_lower_prop:.3f}, {ci_upper_prop:.3f}]")
print(f"95% CI (percentage): [{ci_lower_prop*100:.1f}%, {ci_upper_prop*100:.1f}%]")

# Effect of confidence level
print("\n" + "="*50)
print("EFFECT OF CONFIDENCE LEVEL")
print("="*50)

confidence_levels = [0.90, 0.95, 0.99]
for conf_level in confidence_levels:
alpha = 1 - conf_level
t_crit = stats.t.ppf(1 - alpha/2, df=len(sample_data)-1)
me = t_crit * se
ci_l = sample_mean - me
ci_u = sample_mean + me
print(f"{conf_level*100:.0f}% CI: [{ci_l:.2f}, {ci_u:.2f}], Width: {ci_u-ci_l:.2f}")

print("\nObservation: Higher confidence level → Wider interval")

# Visualization
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Plot 1: Sample means and CI
np.random.seed(42)
sample_means = []
cis_lower = []
cis_upper = []

true_mean = 100
for i in range(50):
sample = np.random.normal(true_mean, 15, 30)
m = sample.mean()
s = sample.std(ddof=1)
se_i = s / np.sqrt(30)
t_crit = stats.t.ppf(0.975, df=29)

sample_means.append(m)
cis_lower.append(m - t_crit * se_i)
cis_upper.append(m + t_crit * se_i)

colors = ['red' if true_mean < cis_lower[i] or true_mean > cis_upper[i] else 'blue'
for i in range(50)]

for i in range(50):
axes[0].plot([cis_lower[i], cis_upper[i]], [i, i], color=colors[i], linewidth=2)
axes[0].scatter(sample_means[i], i, color=colors[i], s=50, zorder=5)

axes[0].axvline(true_mean, color='green', linestyle='--', linewidth=2, label='True Mean')
axes[0].set_xlabel('Value')
axes[0].set_ylabel('Sample Number')
axes[0].set_title('50 Confidence Intervals (95%)
(Blue include true mean, Red do not)')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Plot 2: Effect of sample size
sample_sizes = [10, 30, 50, 100, 200]
margin_errors = []

for n in sample_sizes:
se_n = 15 / np.sqrt(n)
t_crit = stats.t.ppf(0.975, df=n-1)
me = t_crit * se_n
margin_errors.append(me)

axes[1].plot(sample_sizes, margin_errors, marker='o', linewidth=2, markersize=8, color='blue')
axes[1].set_xlabel('Sample Size')
axes[1].set_ylabel('Margin of Error')
axes[1].set_title('Effect of Sample Size on Margin of Error
(95% CI)')
axes[1].grid(True, alpha=0.3)
axes[1].set_xscale('log')

plt.tight_layout()
plt.show()

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

print("="*50)
print("BOOTSTRAPPING")
print("="*50)

print("""
WHAT IS BOOTSTRAPPING?
- Resampling method: repeatedly sample from data with replacement
- Estimates distribution of a statistic without assumptions
- Very powerful for complex statistics where theoretical distributions unknown

STEPS:
1. Draw random sample WITH REPLACEMENT from original sample
2. Calculate statistic of interest on resampled data
3. Repeat steps 1-2 many times (1000-10000)
4. Use distribution of results to estimate CI and SE
""")

# Original data
original_data = np.array([85, 92, 78, 88, 95, 82, 90, 87, 93, 80, 86, 91])
print(f"\nOriginal data: {original_data}")
print(f"Original mean: {original_data.mean():.2f}")
print(f"Original median: {np.median(original_data):.2f}")

# Bootstrap for mean
print("\n" + "="*50)
print("BOOTSTRAP: MEAN")
print("="*50)

n_bootstrap = 10000
bootstrap_means = []

np.random.seed(42)
for i in range(n_bootstrap):
# Resample with replacement
bootstrap_sample = np.random.choice(original_data, len(original_data), replace=True)
bootstrap_means.append(bootstrap_sample.mean())

bootstrap_means = np.array(bootstrap_means)

print(f"\nBootstrap distribution statistics:")
print(f"Mean of bootstrap means: {bootstrap_means.mean():.2f}")
print(f"Std of bootstrap means (SE): {bootstrap_means.std():.2f}")

# 95% CI using percentile method
ci_lower = np.percentile(bootstrap_means, 2.5)
ci_upper = np.percentile(bootstrap_means, 97.5)
print(f"95% CI (percentile method): [{ci_lower:.2f}, {ci_upper:.2f}]")

# Bootstrap for median
print("\n" + "="*50)
print("BOOTSTRAP: MEDIAN")
print("="*50)

bootstrap_medians = []
for i in range(n_bootstrap):
bootstrap_sample = np.random.choice(original_data, len(original_data), replace=True)
bootstrap_medians.append(np.median(bootstrap_sample))

bootstrap_medians = np.array(bootstrap_medians)

print(f"\nBootstrap median statistics:")
print(f"Mean of bootstrap medians: {bootstrap_medians.mean():.2f}")
print(f"Std of bootstrap medians (SE): {bootstrap_medians.std():.2f}")

ci_lower_median = np.percentile(bootstrap_medians, 2.5)
ci_upper_median = np.percentile(bootstrap_medians, 97.5)
print(f"95% CI for median: [{ci_lower_median:.2f}, {ci_upper_median:.2f}]")

# Bootstrap for correlation
print("\n" + "="*50)
print("BOOTSTRAP: CORRELATION")
print("="*50)

# Create paired data
np.random.seed(42)
x = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
y = 2 * x + np.random.normal(0, 2, 10)

original_corr = np.corrcoef(x, y)[0, 1]
print(f"Original correlation: {original_corr:.3f}")

bootstrap_corrs = []
for i in range(n_bootstrap):
indices = np.random.choice(len(x), len(x), replace=True)
x_boot = x[indices]
y_boot = y[indices]
corr_boot = np.corrcoef(x_boot, y_boot)[0, 1]
bootstrap_corrs.append(corr_boot)

bootstrap_corrs = np.array(bootstrap_corrs)

print(f"Bootstrap correlation mean: {bootstrap_corrs.mean():.3f}")
print(f"Bootstrap correlation SE: {bootstrap_corrs.std():.3f}")

ci_lower_corr = np.percentile(bootstrap_corrs, 2.5)
ci_upper_corr = np.percentile(bootstrap_corrs, 97.5)
print(f"95% CI for correlation: [{ci_lower_corr:.3f}, {ci_upper_corr:.3f}]")

# Visualization
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

# Bootstrap means
axes[0].hist(bootstrap_means, bins=50, alpha=0.7, color='blue', edgecolor='black')
axes[0].axvline(original_data.mean(), color='red', linestyle='--', linewidth=2, label='Original Mean')
axes[0].axvline(ci_lower, color='green', linestyle='--', linewidth=2, label='95% CI')
axes[0].axvline(ci_upper, color='green', linestyle='--', linewidth=2)
axes[0].set_title('Bootstrap Distribution of Mean')
axes[0].set_xlabel('Sample Mean')
axes[0].set_ylabel('Frequency')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Bootstrap medians
axes[1].hist(bootstrap_medians, bins=50, alpha=0.7, color='green', edgecolor='black')
axes[1].axvline(np.median(original_data), color='red', linestyle='--', linewidth=2, label='Original Median')
axes[1].axvline(ci_lower_median, color='orange', linestyle='--', linewidth=2, label='95% CI')
axes[1].axvline(ci_upper_median, color='orange', linestyle='--', linewidth=2)
axes[1].set_title('Bootstrap Distribution of Median')
axes[1].set_xlabel('Sample Median')
axes[1].set_ylabel('Frequency')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

# Bootstrap correlations
axes[2].hist(bootstrap_corrs, bins=50, alpha=0.7, color='purple', edgecolor='black')
axes[2].axvline(original_corr, color='red', linestyle='--', linewidth=2, label='Original Correlation')
axes[2].axvline(ci_lower_corr, color='orange', linestyle='--', linewidth=2, label='95% CI')
axes[2].axvline(ci_upper_corr, color='orange', linestyle='--', linewidth=2)
axes[2].set_title('Bootstrap Distribution of Correlation')
axes[2].set_xlabel('Sample Correlation')
axes[2].set_ylabel('Frequency')
axes[2].legend()
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

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

print("="*50)
print("BAYESIAN THINKING")
print("="*50)

print("""
BAYES' THEOREM:
P(A|B) = P(B|A) × P(A) / P(B)

COMPONENTS:
- P(A|B) = Posterior (what we want to know)
- P(B|A) = Likelihood (probability of evidence given hypothesis)
- P(A) = Prior (what we believed before data)
- P(B) = Marginal likelihood (normalizing factor)

KEY INSIGHT:
Bayesian approach updates beliefs as we get new data
""")

# Example 1: Disease testing revisited
print("\nEXAMPLE 1: MEDICAL TESTING")
print("="*50)

# Prior: 1% of population has disease
prior_disease = 0.01
prior_no_disease = 1 - prior_disease

# Test accuracy
sensitivity = 0.95 # P(+|disease)
specificity = 0.90 # P(-|no disease)

# What we observe: positive test
# P(+|disease) × P(disease) = likelihood × prior
prob_positive_given_disease = sensitivity * prior_disease
prob_positive_given_no_disease = (1 - specificity) * prior_no_disease

# P(+) = total probability of positive
prob_positive_total = prob_positive_given_disease + prob_positive_given_no_disease

# Posterior: P(disease|+)
posterior_disease = prob_positive_given_disease / prob_positive_total

print(f"Prior P(disease) = {prior_disease:.3f}")
print(f"Sensitivity P(+|disease) = {sensitivity:.3f}")
print(f"Specificity P(-|no disease) = {specificity:.3f}")
print(f"False positive rate = {1-specificity:.3f}")
print(f"\nAfter positive test:")
print(f"Posterior P(disease|+) = {posterior_disease:.3f} or {posterior_disease*100:.1f}%")
print(f"\nInterpretation: Even with positive test, only {posterior_disease*100:.1f}% chance of disease!")

# Example 2: Coin fairness
print("\n" + "="*50)
print("EXAMPLE 2: IS COIN FAIR?")
print("="*50)

# Prior: assume coin is fair (50-50)
# Data: flip 10 times, get 8 heads

print(f"Prior belief: coin is fair (P(heads)=0.5)")
print(f"Data: 10 flips, 8 heads")

# Calculate likelihood for different values of p
p_values = np.linspace(0, 1, 100)
data = 8 # heads
trials = 10

# Likelihood (binomial)
likelihood = stats.binom.pmf(data, trials, p_values)

# Uniform prior
prior = np.ones_like(p_values)

# Posterior (proportional to likelihood × prior)
posterior = likelihood * prior
posterior = posterior / posterior.sum() # Normalize

# Plot
plt.figure(figsize=(12, 6))
plt.plot(p_values, likelihood, label='Likelihood (Data: 8H in 10 flips)', linewidth=2)
plt.plot(p_values, prior, label='Prior (uniform)', linewidth=2)
plt.plot(p_values, posterior, label='Posterior', linewidth=3)

plt.axvline(0.5, color='red', linestyle='--', label='Fair coin (p=0.5)')
plt.xlabel('Probability of Heads (p)')
plt.ylabel('Probability Density')
plt.title('Bayesian Inference: Is the Coin Fair?')
plt.legend(fontsize=11)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# Calculate posterior mean
posterior_mean = np.average(p_values, weights=posterior)
print(f"\nPosterior mean (best estimate): {posterior_mean:.3f}")
print(f"Our estimate shifted from 0.5 (prior) to {posterior_mean:.3f} (posterior)")

# Example 3: Updating beliefs
print("\n" + "="*50)
print("EXAMPLE 3: SEQUENTIAL LEARNING")
print("="*50)

print(f"Starting belief: coin is fair (p=0.5)")
print(f"\nSequence of events:")

p = 0.8 # True probability
prior_mean = 0.5

observations = [1, 1, 0, 1, 1, 1, 0, 1, 1, 1] # 1=heads, 0=tails

for i, obs in enumerate(observations, 1):
# Count heads so far
heads_count = sum(observations[:i])

# Update using Beta-Binomial model
# Beta distribution for prior on p
alpha = 1 + heads_count
beta = 1 + (i - heads_count)

# Posterior mean
posterior_mean_beta = alpha / (alpha + beta)

outcome = "H" if obs == 1 else "T"
print(f"Flip {i:2d}: {outcome} → Heads: {heads_count}, Posterior mean: {posterior_mean_beta:.3f}")

print(f"\nFinal posterior mean: {posterior_mean_beta:.3f}")
print(f"(True value was p={p})")

# Compare Frequentist vs Bayesian
print("\n" + "="*50)
print("FREQUENTIST vs BAYESIAN")
print("="*50)

print(f"""
FREQUENTIST:
- Probability is long-run frequency
- Parameter is fixed, unknown
- Uses only the data
- P-values and confidence intervals

BAYESIAN:
- Probability is subjective degree of belief
- Parameter has a distribution
- Combines data with prior knowledge
- Posterior distributions and credible intervals

WHEN TO USE:
- Frequentist: When you have lots of data, want standard tests
- Bayesian: When you want to incorporate prior knowledge, small samples
""")