统计推断

import numpy as np

# 从正态总体中抽样
np.random.seed(42)
population_mean = 100
population_std = 15
sample_size = 30

sample = np.random.normal(population_mean, population_std, sample_size)
sample_mean = np.mean(sample)
sample_var = np.var(sample, ddof=1)  # Bessel 校正
print(f"样本均值: {sample_mean:.2f}, 样本方差: {sample_var:.2f}")

import scipy.stats as stats

# 验证中心极限定理
n_simulations = 10000
sample_means = []
for _ in range(n_simulations):
    s = np.random.exponential(scale=2, size=50)
    sample_means.append(np.mean(s))

# 样本均值的分布应近似正态
mean_of_means = np.mean(sample_means)
std_of_means = np.std(sample_means)
print(f"CLT 验证: 均值={mean_of_means:.4f}, 标准差={std_of_means:.4f}")

import numpy as np

# 矩估计 vs 样本统计量
data = np.array([12.5, 14.2, 11.8, 13.6, 15.1, 12.9, 14.7, 13.3])

# 一阶矩估计 = 样本均值
theta_hat_mean = np.mean(data)
# 二阶矩估计方差
theta_hat_var = np.mean((data - theta_hat_mean) ** 2)
# 无偏方差（Bessel 校正）
theta_hat_var_unbiased = np.var(data, ddof=1)

print(f"矩估计均值: {theta_hat_mean:.3f}")
print(f"矩估计方差: {theta_hat_var:.3f}")
print(f"无偏方差:   {theta_hat_var_unbiased:.3f}")

from scipy import stats

# 评估估计量的无偏性和有效性
def simulate_estimator(n, true_theta, n_sims=5000):
    estimates = []
    for _ in range(n_sims):
        sample = np.random.exponential(scale=1/true_theta, size=n)
        estimates.append(1 / np.mean(sample))  # 指数分布λ的MLE
    return {
        "bias": np.mean(estimates) - true_theta,
        "variance": np.var(estimates),
        "mse": np.mean((np.array(estimates) - true_theta)**2)
    }

result = simulate_estimator(n=20, true_theta=2.0)
print(f"偏差: {result['bias']:.4f}, 方差: {result['variance']:.4f}, MSE: {result['mse']:.4f}")

import numpy as np
from scipy.optimize import minimize

# 正态分布 MLE（解析解 vs 数值解）
data = np.random.normal(loc=5, scale=2, size=200)

# 解析解：MLE 均值 = 样本均值, MLE 方差 = 样本方差(n 分母)
mu_mle = np.mean(data)
sigma2_mle = np.mean((data - mu_mle)**2)
print(f"解析解: μ={mu_mle:.3f}, σ²={sigma2_mle:.3f}")

# 数值优化求解
def neg_log_likelihood(params):
    mu, sigma = params
    if sigma <= 0:
        return 1e10
    n = len(data)
    return n/2 * np.log(2*np.pi*sigma2) + np.sum((data-mu)2)/(2*sigma**2)

result = minimize(neg_log_likelihood, x0=[0, 1], method="Nelder-Mead")
print(f"数值解: μ={result.x[0]:.3f}, σ={result.x[1]:.3f}")

# 伯努利分布的 MLE：估计硬币正面概率
def bernoulli_mle(data):
    """二项分布中成功概率 p 的 MLE = 样本均值"""
    return np.mean(data)

# 模拟：抛 100 次硬币，正面 62 次
coin_flips = np.concatenate([np.ones(62), np.zeros(38)])
p_hat = bernoulli_mle(coin_flips)
print(f"MLE 估计正面概率: {p_hat:.3f}")

# Fisher 信息量
I_n = len(coin_flips) / (p_hat * (1 - p_hat))
print(f"Fisher 信息量: {I_n:.2f}")
print(f"渐近方差: {1/I_n:.6f}")

import numpy as np
from scipy.stats import norm

# 正态先验下的 MAP 估计
# 数据: X ~ N(θ, σ²), 先验: θ ~ N(μ₀, τ²)
data = np.array([3.2, 4.1, 2.8, 3.9, 4.5, 3.1, 3.7, 4.3])
sigma2 = 1.0   # 已知数据方差
mu0 = 0.0      # 先验均值
tau2 = 10.0    # 先验方差（大=弱先验）

# MAP 解析解: θ̂ = (n·X̄/σ² + μ₀/τ²) / (n/σ² + 1/τ²)
n = len(data)
x_bar = np.mean(data)
theta_map = (n * x_bar / sigma2 + mu0 / tau2) / (n / sigma2 + 1 / tau2)
theta_mle = x_bar

print(f"MLE: {theta_mle:.3f}")
print(f"MAP: {theta_map:.3f}")
print(f"先验越强(MAP越靠近μ₀)，tau2 越小效果越明显")

# MAP 与正则化的等价性：Logistic 回归
from sklearn.linear_model import LogisticRegression
import numpy as np

np.random.seed(42)
X = np.random.randn(200, 5)
w_true = np.array([2, -1, 0.5, 0, 0])
y = (X @ w_true + np.random.randn(200) * 0.5 > 0).astype(int)

# L2 正则化 = 高斯先验的 MAP
model_l2 = LogisticRegression(penalty="l2", C=1.0, solver="lbfgs")
model_l2.fit(X, y)

# 无正则化 = MLE
model_mle = LogisticRegression(penalty=None, solver="lbfgs")
model_mle.fit(X, y)

print("MAP(L2) 系数:", np.round(model_l2.coef_[0], 3))
print("MLE 系数:   ", np.round(model_mle.coef_[0], 3))

import numpy as np
import scipy.stats as stats

# t 分布置信区间（小样本，σ 未知）
data = np.array([23.1, 25.4, 22.8, 24.6, 26.1, 23.9, 25.2, 24.3, 23.7, 25.8])
n = len(data)
x_bar = np.mean(data)
s = np.std(data, ddof=1)
alpha = 0.05

# t 临界值
t_crit = stats.t.ppf(1 - alpha/2, df=n-1)
margin = t_crit * s / np.sqrt(n)
ci_lower = x_bar - margin
ci_upper = x_bar + margin

print(f"95% CI: [{ci_lower:.3f}, {ci_upper:.3f}]")
print(f"点估计: {x_bar:.3f} ± {margin:.3f}")

# Bootstrap 置信区间（无分布假设）
def bootstrap_ci(data, n_bootstrap=10000, alpha=0.05):
    """Bootstrap 百分位数法计算置信区间"""
    boot_means = []
    for _ in range(n_bootstrap):
        resample = np.random.choice(data, size=len(data), replace=True)
        boot_means.append(np.mean(resample))
    
    lower = np.percentile(boot_means, 100 * alpha / 2)
    upper = np.percentile(boot_means, 100 * (1 - alpha / 2))
    return lower, upper

data = np.array([12, 15, 8, 22, 18, 11, 25, 14, 19, 16, 13, 21])
ci = bootstrap_ci(data)
print(f"Bootstrap 95% CI: [{ci[0]:.3f}, {ci[1]:.3f}]")
print(f"传统 95% CI: 请对比与 t 分布方法的差异")

import numpy as np
import scipy.stats as stats

# 单样本 t 检验：新药平均效果是否 ≠ 0
np.random.seed(42)
drug_effects = np.array([2.3, -0.5, 1.8, 3.1, 0.9, -1.2, 2.7, 1.5, 0.3, 2.9])

t_stat, p_value = stats.ttest_1samp(drug_effects, popmean=0)
print(f"t 统计量: {t_stat:.4f}")
print(f"p 值: {p_value:.4f}")
print(f"α=0.05 下: {'拒绝 H₀' if p_value < 0.05 else '不拒绝 H₀'}")

# 效应量 (Cohen's d)
cohens_d = np.mean(drug_effects) / np.std(drug_effects, ddof=1)
print(f"Cohen's d: {cohens_d:.3f}")

# 功效分析：需要多少样本？
from statsmodels.stats.power import TTestIndPower

analysis = TTestIndPower()

# 计算功效（已知效应量和样本量）
effect_size = 0.5  # 中等效应
n_per_group = 64
power = analysis.power(effect_size=effect_size, nobs1=n_per_group, alpha=0.05)
print(f"效应量={effect_size}, 每组 n={n_per_group}, 功效={power:.3f}")

# 计算所需样本量（已知期望功效）
n_needed = analysis.solve_power(effect_size=effect_size, power=0.8, alpha=0.05)
print(f"达到 80% 功效需要每组 n={int(np.ceil(n_needed))}")

import numpy as np
from scipy import stats

# 贝叶斯 A/B 测试：Beta-Binomial 共轭
# 版本 A: 1200 次展示, 96 次转化; 版本 B: 1350 次展示, 128 次转化
alpha_prior, beta_prior = 1, 1  # 均匀先验 Beta(1,1)

# 后验参数
alpha_a = alpha_prior + 96
beta_a = beta_prior + 1200 - 96
alpha_b = alpha_prior + 128
beta_b = beta_prior + 1350 - 128

# Monte Carlo 估计 P(B > A)
samples_a = np.random.beta(alpha_a, beta_a, 100000)
samples_b = np.random.beta(alpha_b, beta_b, 100000)
prob_b_better = np.mean(samples_b > samples_a)
lift = np.mean(samples_b - samples_a)

print(f"P(B > A | 数据) = {prob_b_better:.4f}")
print(f"转化率提升: {lift:.4f} ({lift/np.mean(samples_a)*100:.1f}%)")
print(f"A 的后验均值: {stats.beta.mean(alpha_a, beta_a):.4f}")
print(f"B 的后验均值: {stats.beta.mean(alpha_b, beta_b):.4f}")

# 贝叶斯线性回归
np.random.seed(42)
n_train = 30
X = np.linspace(0, 10, n_train).reshape(-1, 1)
X_design = np.column_stack([np.ones(n_train), X])
y = 2 * X.flatten() + 1 + np.random.randn(n_train) * 2

# 先验: w ~ N(0, σ_w²·I), 噪声: y|X,w ~ N(Xw, σ²)
sigma2 = 4.0
sigma_w2 = 10.0
w_prior_mean = np.zeros(2)
w_prior_cov = sigma_w2 * np.eye(2)

# 后验: w|X,y ~ N(w_n, Σ_n)
XtX = X_design.T @ X_design
w_post_cov = np.linalg.inv(np.linalg.inv(w_prior_cov) + XtX / sigma2)
w_post_mean = w_post_cov @ (
    np.linalg.inv(w_prior_cov) @ w_prior_mean + X_design.T @ y / sigma2
)

print(f"后验权重: w₀={w_post_mean[0]:.3f}, w₁={w_post_mean[1]:.3f}")

# 预测不确定性
x_new = np.array([[1, 5.5]])
pred_mean = x_new @ w_post_mean
pred_var = x_new @ w_post_cov @ x_new.T + sigma2
print(f"x=5.5 预测: {pred_mean[0]:.3f} ± {1.96*np.sqrt(pred_var[0,0]):.3f}")