Statistics Interview Questions - Hard

 1import numpy as np
 2from scipy import stats
 3from statsmodels.stats.multitest import multipletests
 4
 5np.random.seed(42)
 6
 7# Simulate multiple hypothesis tests
 8n_tests = 20
 9alpha = 0.05
10
11# Generate p-values (most from null, few from alternative)
12p_values = []
13true_hypotheses = []  # Track which are actually true
14
15for i in range(n_tests):
16    if i < 16:  # 16 true null hypotheses
17        # Null hypothesis: mean = 0
18        data = np.random.normal(0, 1, 30)
19        true_hypotheses.append(True)
20    else:  # 4 false null hypotheses (true effects exist)
21        # Alternative: mean = 0.8
22        data = np.random.normal(0.8, 1, 30)
23        true_hypotheses.append(False)
24    
25    _, p_val = stats.ttest_1samp(data, 0)
26    p_values.append(p_val)
27
28p_values = np.array(p_values)
29true_hypotheses = np.array(true_hypotheses)
30
31print("Multiple Comparisons Problem:")
32print(f"Number of tests: {n_tests}")
33print(f"Nominal α: {alpha}")
34print(f"Expected false positives without correction: {alpha * n_tests:.1f}")
35
36# Uncorrected
37rejections_uncorrected = p_values < alpha
38false_positives_uncorrected = np.sum(rejections_uncorrected & true_hypotheses)
39print(f"\nUncorrected:")
40print(f"  Rejections: {np.sum(rejections_uncorrected)}")
41print(f"  False positives: {false_positives_uncorrected}")
42
43# Bonferroni correction
44alpha_bonf = alpha / n_tests
45rejections_bonf = p_values < alpha_bonf
46false_positives_bonf = np.sum(rejections_bonf & true_hypotheses)
47print(f"\nBonferroni correction (α = {alpha_bonf:.4f}):")
48print(f"  Rejections: {np.sum(rejections_bonf)}")
49print(f"  False positives: {false_positives_bonf}")
50print(f"  Trade-off: Very conservative, may miss real effects")
51
52# FDR (Benjamini-Hochberg)
53rejections_fdr, p_corrected_fdr, alpha_sidak, alpha_bonf_fdr = multipletests(
54    p_values, alpha=alpha, method='fdr_bh'
55)
56false_positives_fdr = np.sum(rejections_fdr & true_hypotheses)
57print(f"\nFalse Discovery Rate (FDR) - Benjamini-Hochberg:")
58print(f"  Rejections: {np.sum(rejections_fdr)}")
59print(f"  False positives: {false_positives_fdr}")
60print(f"  Trade-off: Less conservative, controls FDR instead of FWER")
61
62# Holm-Bonferroni (step-down)
63rejections_holm, p_corrected_holm, _, _ = multipletests(
64    p_values, alpha=alpha, method='holm'
65)
66false_positives_holm = np.sum(rejections_holm & true_hypotheses)
67print(f"\nHolm-Bonferroni (step-down):")
68print(f"  Rejections: {np.sum(rejections_holm)}")
69print(f"  False positives: {false_positives_holm}")
70print(f"  Trade-off: Less conservative than Bonferroni, still controls FWER")
71
72# Comparison
73print(f"\nSummary:")
74print(f"{'Method':<20} {'Rejections':<12} {'False Positives':<15}")
75print("-" * 50)
76print(f"{'Uncorrected':<20} {np.sum(rejections_uncorrected):<12} {false_positives_uncorrected:<15}")
77print(f"{'Bonferroni':<20} {np.sum(rejections_bonf):<12} {false_positives_bonf:<15}")
78print(f"{'FDR (BH)':<20} {np.sum(rejections_fdr):<12} {false_positives_fdr:<15}")
79print(f"{'Holm':<20} {np.sum(rejections_holm):<12} {false_positives_holm:<15}")

 1import numpy as np
 2import matplotlib.pyplot as plt
 3from scipy import stats
 4from scipy.optimize import minimize
 5
 6# Example: MLE for normal distribution parameters
 7np.random.seed(42)
 8
 9# True parameters
10true_mu = 5.0
11true_sigma = 2.0
12n = 100
13
14# Generate data
15data = np.random.normal(true_mu, true_sigma, n)
16
17# MLE for normal distribution (analytical solution)
18# μ_MLE = sample mean
19# σ²_MLE = sample variance (biased version)
20mu_mle = np.mean(data)
21sigma_mle = np.std(data, ddof=0)  # MLE uses N, not N-1
22
23print("MLE for Normal Distribution:")
24print(f"True μ = {true_mu:.3f}, True σ = {true_sigma:.3f}")
25print(f"MLE μ = {mu_mle:.3f}, MLE σ = {sigma_mle:.3f}")
26
27# Numerical MLE using optimization
28def negative_log_likelihood(params, data):
29    """Negative log-likelihood for normal distribution."""
30    mu, sigma = params
31    if sigma <= 0:
32        return np.inf
33    return -np.sum(stats.norm.logpdf(data, loc=mu, scale=sigma))
34
35# Optimize
36initial_guess = [0, 1]
37result = minimize(negative_log_likelihood, initial_guess, args=(data,), 
38                  method='L-BFGS-B', bounds=[(None, None), (0.01, None)])
39
40mu_mle_opt = result.x[0]
41sigma_mle_opt = result.x[1]
42
43print(f"\nNumerical MLE:")
44print(f"μ = {mu_mle_opt:.3f}, σ = {sigma_mle_opt:.3f}")
45print(f"Log-likelihood: {-result.fun:.3f}")
46
47# MLE for Poisson distribution
48true_lambda = 3.0
49poisson_data = np.random.poisson(true_lambda, n)
50
51# Analytical MLE for Poisson: λ_MLE = sample mean
52lambda_mle = np.mean(poisson_data)
53
54print(f"\nMLE for Poisson Distribution:")
55print(f"True λ = {true_lambda:.3f}")
56print(f"MLE λ = {lambda_mle:.3f}")
57
58# Likelihood function visualization
59def log_likelihood_poisson(lambda_val, data):
60    """Log-likelihood for Poisson."""
61    return np.sum(stats.poisson.logpmf(data, lambda_val))
62
63lambda_vals = np.linspace(1, 5, 100)
64log_likes = [log_likelihood_poisson(l, poisson_data) for l in lambda_vals]
65
66# Find MLE
67mle_idx = np.argmax(log_likes)
68lambda_mle_vis = lambda_vals[mle_idx]
69
70print(f"\nMLE from likelihood plot: {lambda_mle_vis:.3f}")
71
72# Asymptotic distribution of MLE
73# For large n, MLE ~ Normal(θ, I(θ)⁻¹)
74# For Poisson: I(λ) = n/λ, so var(λ_MLE) ≈ λ/n
75n_samples = 1000
76mle_samples = [np.mean(np.random.poisson(true_lambda, n)) for _ in range(n_samples)]
77
78theoretical_var = true_lambda / n
79theoretical_std = np.sqrt(theoretical_var)
80empirical_std = np.std(mle_samples)
81
82print(f"\nAsymptotic Distribution of MLE (Poisson):")
83print(f"Theoretical std: {theoretical_std:.4f}")
84print(f"Empirical std: {empirical_std:.4f}")
85print(f"Sample mean (should be close to true λ): {np.mean(mle_samples):.3f}")
86
87# Verify asymptotic normality
88from scipy.stats import normaltest
89stat, p_val = normaltest(mle_samples)
90print(f"Normality test p-value: {p_val:.4f}")

 1import numpy as np
 2from scipy import stats
 3import matplotlib.pyplot as plt
 4
 5np.random.seed(42)
 6
 7# Original data (skewed, non-normal)
 8data = np.random.exponential(scale=2, size=30)
 9true_median = 2 * np.log(2)  # True median of exponential
10
11print("Bootstrap Confidence Interval Example:")
12print(f"Sample size: {len(data)}")
13print(f"Sample median: {np.median(data):.3f}")
14print(f"True median: {true_median:.3f}")
15
16# Bootstrap procedure
17n_bootstrap = 10000
18bootstrap_medians = []
19
20for _ in range(n_bootstrap):
21    # Resample with replacement
22    bootstrap_sample = np.random.choice(data, size=len(data), replace=True)
23    bootstrap_medians.append(np.median(bootstrap_sample))
24
25bootstrap_medians = np.array(bootstrap_medians)
26
27# Confidence interval using percentile method
28ci_lower = np.percentile(bootstrap_medians, 2.5)
29ci_upper = np.percentile(bootstrap_medians, 97.5)
30
31print(f"\nBootstrap 95% CI (percentile method): [{ci_lower:.3f}, {ci_upper:.3f}]")
32
33# Bootstrap standard error
34bootstrap_se = np.std(bootstrap_medians)
35print(f"Bootstrap standard error: {bootstrap_se:.3f}")
36
37# Compare with normal approximation (may not be valid)
38sample_median = np.median(data)
39# Bootstrap t-interval
40t_critical = stats.t.ppf(0.975, df=len(data)-1)
41ci_normal = [sample_median - t_critical * bootstrap_se, 
42             sample_median + t_critical * bootstrap_se]
43print(f"Bootstrap 95% CI (t-interval): [{ci_normal[0]:.3f}, {ci_normal[1]:.3f}]")
44
45# Bootstrap for correlation (no simple CI formula)
46np.random.seed(42)
47x = np.random.normal(0, 1, 50)
48y = 2 + 1.5*x + np.random.normal(0, 0.5, 50)
49true_correlation = 1.5 / np.sqrt(1 + 0.5**2)  # Theoretical correlation
50
51sample_correlation = np.corrcoef(x, y)[0, 1]
52print(f"\nBootstrap for Correlation:")
53print(f"Sample correlation: {sample_correlation:.3f}")
54print(f"Theoretical correlation: {true_correlation:.3f}")
55
56bootstrap_correlations = []
57for _ in range(n_bootstrap):
58    indices = np.random.choice(len(x), size=len(x), replace=True)
59    boot_x = x[indices]
60    boot_y = y[indices]
61    boot_corr = np.corrcoef(boot_x, boot_y)[0, 1]
62    bootstrap_correlations.append(boot_corr)
63
64bootstrap_correlations = np.array(bootstrap_correlations)
65ci_corr = np.percentile(bootstrap_correlations, [2.5, 97.5])
66print(f"Bootstrap 95% CI for correlation: [{ci_corr[0]:.3f}, {ci_corr[1]:.3f}]")
67
68# Bootstrap hypothesis testing
69# Test: H0: median = 1.0
70null_value = 1.0
71test_statistic = np.median(data) - null_value
72
73# Bootstrap under H0: center data at null value
74centered_data = data - np.median(data) + null_value
75bootstrap_test_stats = []
76
77for _ in range(n_bootstrap):
78    boot_sample = np.random.choice(centered_data, size=len(data), replace=True)
79    boot_stat = np.median(boot_sample) - null_value
80    bootstrap_test_stats.append(boot_stat)
81
82bootstrap_test_stats = np.array(bootstrap_test_stats)
83
84# Two-tailed p-value
85p_value = np.mean(np.abs(bootstrap_test_stats) >= np.abs(test_statistic))
86print(f"\nBootstrap Hypothesis Test:")
87print(f"H0: median = {null_value}")
88print(f"Test statistic: {test_statistic:.3f}")
89print(f"Bootstrap p-value: {p_value:.3f}")

 1import numpy as np
 2from scipy import stats
 3import matplotlib.pyplot as plt
 4
 5# Bayesian inference for normal mean (known variance)
 6
 7# Data
 8np.random.seed(42)
 9true_mu = 10
10true_sigma = 2
11n = 20
12data = np.random.normal(true_mu, true_sigma, n)
13sample_mean = np.mean(data)
14
15print("Bayesian Inference: Normal Mean (known variance)")
16print(f"True μ = {true_mu}, True σ = {true_sigma}")
17print(f"Sample mean = {sample_mean:.3f}, n = {n}")
18
19# Conjugate prior: Normal(μ₀, τ₀²)
20# Prior
21mu_0 = 8  # Prior mean
22tau_0_sq = 4  # Prior variance
23tau_0 = np.sqrt(tau_0_sq)
24
25# Posterior (conjugate prior)
26# Posterior is Normal with:
27# μ_n = (τ₀² * x̄ + σ²/n * μ₀) / (τ₀² + σ²/n)
28# τ_n² = 1 / (1/τ₀² + n/σ²)
29
30tau_n_sq = 1 / (1/tau_0_sq + n / (true_sigma**2))
31mu_n = (tau_0_sq * sample_mean + (true_sigma**2/n) * mu_0) / (tau_0_sq + true_sigma**2/n)
32tau_n = np.sqrt(tau_n_sq)
33
34print(f"\nPrior: Normal(μ₀={mu_0}, τ₀²={tau_0_sq})")
35print(f"Posterior: Normal(μ_n={mu_n:.3f}, τ_n²={tau_n_sq:.3f})")
36
37# Compare with frequentist
38frequentist_mean = sample_mean
39frequentist_se = true_sigma / np.sqrt(n)
40frequentist_ci = [frequentist_mean - 1.96*frequentist_se, 
41                   frequentist_mean + 1.96*frequentist_se]
42
43bayesian_ci = stats.norm.interval(0.95, loc=mu_n, scale=tau_n)
44
45print(f"\nFrequentist 95% CI: [{frequentist_ci[0]:.3f}, {frequentist_ci[1]:.3f}]")
46print(f"Bayesian 95% Credible Interval: [{bayesian_ci[0]:.3f}, {bayesian_ci[1]:.3f}]")
47
48# Visualize prior, likelihood, posterior
49mu_range = np.linspace(5, 15, 200)
50
51prior = stats.norm.pdf(mu_range, loc=mu_0, scale=tau_0)
52# Likelihood (as function of mu)
53likelihood = np.exp(-n/(2*true_sigma**2) * (mu_range - sample_mean)**2)
54likelihood = likelihood / np.trapz(likelihood, mu_range) * np.max(prior) / np.max(likelihood)
55posterior = stats.norm.pdf(mu_range, loc=mu_n, scale=tau_n)
56
57# Bayesian A/B testing (Beta-Binomial)
58print("\n" + "="*50)
59print("Bayesian A/B Testing")
60
61# Priors: Beta(α, β) - uniform priors
62alpha_A, beta_A = 1, 1
63alpha_B, beta_B = 1, 1
64
65# Data
66successes_A, trials_A = 45, 100
67successes_B, trials_B = 55, 100
68
69# Posteriors (Beta is conjugate prior for binomial)
70alpha_A_post = alpha_A + successes_A
71beta_A_post = beta_A + (trials_A - successes_A)
72alpha_B_post = alpha_B + successes_B
73beta_B_post = beta_B + (trials_B - successes_B)
74
75# Probability that B > A
76n_samples = 10000
77samples_A = np.random.beta(alpha_A_post, beta_A_post, n_samples)
78samples_B = np.random.beta(alpha_B_post, beta_B_post, n_samples)
79prob_B_better = np.mean(samples_B > samples_A)
80
81print(f"P(B > A) = {prob_B_better:.3f}")
82
83# Expected values
84expected_A = alpha_A_post / (alpha_A_post + beta_A_post)
85expected_B = alpha_B_post / (alpha_B_post + beta_B_post)
86print(f"Expected conversion A: {expected_A:.3f}")
87print(f"Expected conversion B: {expected_B:.3f}")
88
89# Credible intervals
90ci_A = stats.beta.interval(0.95, alpha_A_post, beta_A_post)
91ci_B = stats.beta.interval(0.95, alpha_B_post, beta_B_post)
92print(f"95% Credible Interval A: [{ci_A[0]:.3f}, {ci_A[1]:.3f}]")
93print(f"95% Credible Interval B: [{ci_B[0]:.3f}, {ci_B[1]:.3f}]")
94
95# Advantage: Direct probability statements
96print(f"\nBayesian advantage:")
97print(f"We can say: 'There is a {prob_B_better:.1%} probability that B is better'")
98print(f"Frequentist can only say: 'We reject/fail to reject H0'")

  1import numpy as np
  2import pandas as pd
  3from scipy import stats
  4from itertools import product
  5
  6# Design optimization: Sample size for different designs
  7
  8def calculate_power_two_sample(effect_size, n_per_group, std_dev=1, alpha=0.05):
  9    """Calculate power for two-sample t-test."""
 10    z_alpha = stats.norm.ppf(1 - alpha/2)
 11    z_beta = stats.norm.ppf(0.80)  # Target power
 12    n_required = 2 * ((z_alpha + z_beta)**2 * std_dev**2) / effect_size**2
 13    return 1 - stats.norm.cdf(z_alpha - effect_size * np.sqrt(n_per_group) / std_dev)
 14
 15print("Experimental Design Optimization:")
 16print("="*50)
 17
 18# Scenario: Testing new drug vs placebo
 19effect_size = 0.5  # Standardized effect size
 20alpha = 0.05
 21target_power = 0.80
 22
 23# Required sample size
 24z_alpha = stats.norm.ppf(1 - alpha/2)
 25z_beta = stats.norm.ppf(target_power)
 26n_per_group = int(np.ceil(2 * ((z_alpha + z_beta)**2) / effect_size**2))
 27
 28print(f"Scenario: Drug vs Placebo")
 29print(f"Effect size (Cohen's d): {effect_size}")
 30print(f"Required sample size per group: {n_per_group}")
 31print(f"Total sample size: {n_per_group * 2}")
 32
 33# Compare designs
 34designs = {
 35    'CRD': {'n': n_per_group * 2, 'efficiency': 1.0},
 36    'RBD (4 blocks)': {'n': n_per_group * 2, 'efficiency': 1.2},  # Blocking reduces variance
 37    'Crossover': {'n': int(n_per_group * 0.6), 'efficiency': 1.5}  # Within-subject more efficient
 38}
 39
 40print(f"\nDesign Comparison:")
 41print(f"{'Design':<20} {'Total N':<12} {'Efficiency':<12}")
 42print("-" * 45)
 43for design, props in designs.items():
 44    print(f"{design:<20} {props['n']:<12} {props['efficiency']:<12}")
 45
 46# Factorial design: Testing 2 factors simultaneously
 47print(f"\nFactorial Design (2×2):")
 48print("Testing: Drug (A/B) × Dose (Low/High)")
 49
 50# Main effects and interactions
 51effects = {
 52    'Main Effect A': 0.3,
 53    'Main Effect Dose': 0.4,
 54    'Interaction': 0.2
 55}
 56
 57# Calculate power for each effect
 58n_per_cell = 25
 59for effect_name, effect in effects.items():
 60    power = calculate_power_two_sample(effect, n_per_cell)
 61    print(f"{effect_name}: effect={effect:.2f}, power={power:.3f}")
 62
 63# Optimal allocation for stratified design
 64print(f"\nOptimal Stratification:")
 65strata_sizes = {'Young': 1000, 'Middle': 2000, 'Old': 1000}
 66total_N = sum(strata_sizes.values())
 67target_n = 300
 68
 69# Proportional allocation
 70proportional = {k: int(v/total_N * target_n) for k, v in strata_sizes.items()}
 71print(f"Proportional allocation: {proportional}")
 72
 73# Optimal allocation (Neyman allocation - accounts for variance)
 74# Assume older group has higher variance
 75strata_vars = {'Young': 1.0, 'Middle': 1.2, 'Old': 1.5}
 76weights = {k: v * np.sqrt(strata_vars[k]) for k, v in strata_sizes.items()}
 77total_weight = sum(weights.values())
 78optimal = {k: int(v/total_weight * target_n) for k, v in weights.items()}
 79print(f"Optimal allocation (Neyman): {optimal}")
 80
 81# Randomized block design simulation
 82print(f"\nRandomized Block Design Simulation:")
 83np.random.seed(42)
 84
 85n_blocks = 4
 86n_per_block = 10
 87block_effects = np.random.normal(0, 2, n_blocks)  # Block-to-block variation
 88treatment_effect = 3
 89
 90# Generate data
 91data = []
 92for block in range(n_blocks):
 93    control = np.random.normal(block_effects[block], 1, n_per_block // 2)
 94    treatment = np.random.normal(block_effects[block] + treatment_effect, 1, n_per_block // 2)
 95    
 96    for val in control:
 97        data.append({'block': block, 'treatment': 'control', 'value': val})
 98    for val in treatment:
 99        data.append({'block': block, 'treatment': 'treatment', 'value': val})
100
101df = pd.DataFrame(data)
102
103# Two approaches: Ignore blocks vs account for blocks
104# Ignore blocks (treat as CRD)
105control_vals = df[df['treatment']=='control']['value']
106treat_vals = df[df['treatment']=='treatment']['value']
107t_stat_crd, p_crd = stats.ttest_ind(treat_vals, control_vals)
108
109# Account for blocks (two-way ANOVA or paired t-test per block)
110from scipy.stats import f_oneway
111# Simplified: block as factor
112block_means_control = df[df['treatment']=='control'].groupby('block')['value'].mean()
113block_means_treat = df[df['treatment']=='treatment'].groupby('block')['value'].mean()
114t_stat_rbd, p_rbd = stats.ttest_rel(block_means_treat, block_means_control)
115
116print(f"CRD analysis (ignoring blocks): t={t_stat_crd:.3f}, p={p_crd:.4f}")
117print(f"RBD analysis (accounting for blocks): t={t_stat_rbd:.3f}, p={p_rbd:.4f}")
118print(f"Blocking improves precision (lower p-value)")