统计计算python demo代码

嗯，用户之前已经询问了总离差平方和以及大学和硕士阶段的统计学概念，现在他们想要用Python实现统计计算的相关中间件的示例代码。我需要先回顾一下之前的对话内容，确保提供的代码示例覆盖用户提到的各个模块，比如参数估计、假设检验、回归分析、多元统计分析等。

首先，用户提到的中间件可能指的是Python中的相关库，比如statsmodels、scikit-learn、numpy、pandas等。需要确保每个示例都使用这些库，并且涵盖不同的统计方法。例如，参数估计部分可以使用statsmodels进行线性回归，假设检验可以用scipy.stats中的t检验，非参数检验可能用到Mann-Whitney U检验等。

接下来，要确定用户的需求层次。用户可能是学生或刚入门的数据分析者，希望通过示例代码学习如何应用统计学概念。因此，代码需要简洁明了，注释清晰，并且包含数据生成和结果解释的部分。例如，在展示线性回归时，可以生成模拟数据，然后拟合模型并输出摘要，这样用户可以看到整个过程。

还要注意覆盖用户提到的各个模块，比如时间序列中的ARIMA，生存分析中的Kaplan-Meier曲线，以及机器学习中的随机森林等。每个部分需要独立的代码块，并确保库的导入和版本兼容性。例如，使用pmdarima进行ARIMA模型选择，而lifelines处理生存分析。

另外，用户可能对代码的实际应用场景感兴趣，所以在每个示例前简要说明应用场景，帮助用户理解何时使用该方法。比如，在展示Bootstrap方法时，说明其用于置信区间的估计，适用于小样本数据。

需要检查是否有遗漏的中间件或方法。例如，贝叶斯推断部分可以使用PyMC3，但用户可能没有提到，但根据之前的回答，PyMC3被列在统计计算工具中，所以应该包含一个贝叶斯线性回归的示例。

此外，代码的可运行性很重要。确保所有示例代码都是自包含的，包括必要的库安装（如注释中的安装命令），但实际代码中不需要!pip install，除非用户指定。可能需要提醒用户安装相关库，比如pmdarima和lifelines。

最后，总结部分要强调练习和探索的重要性，并鼓励用户修改参数和数据集来加深理解。同时，提供进一步学习的建议，比如官方文档和实际项目应用。

---

以下是用 Python 实现统计学核心方法的代码示例，涵盖参数估计、假设检验、回归分析、时间序列等场景，使用 statsmodels、scikit-learn 等常用库：

---

1. 参数估计与假设检验

场景：计算均值的置信区间 & 独立样本 t 检验

import numpy as np
import scipy.stats as stats
import statsmodels.stats.api as sms

# 生成模拟数据
np.random.seed(42)
data_A = np.random.normal(5.0, 1.5, 30)
data_B = np.random.normal(5.5, 1.5, 30)

# --- 计算均值的 95% 置信区间 ---
confidence_level = 0.95
ci = stats.t.interval(confidence_level, len(data_A)-1, 
                     loc=np.mean(data_A), 
                     scale=stats.sem(data_A))
print(f"均值置信区间: {ci}")

# --- 独立样本 t 检验 ---
t_stat, p_value = stats.ttest_ind(data_A, data_B)
print(f"t 统计量: {t_stat:.3f}, p 值: {p_value:.3f}")

---

2. 线性回归分析

场景：简单线性回归与诊断

import statsmodels.api as sm
import pandas as pd

# 生成数据
X = np.random.rand(100, 1) * 10  # 特征
y = 2.5 * X.squeeze() + np.random.normal(0, 2, 100)  # 目标变量 + 噪声

# 添加截距项并拟合模型
X = sm.add_constant(X)
model = sm.OLS(y, X).fit()

# 输出回归结果
print(model.summary())

# 残差诊断（可视化）
import matplotlib.pyplot as plt
residuals = model.resid
fig, ax = plt.subplots(1, 2, figsize=(12, 4))
ax[0].scatter(model.fittedvalues, residuals)
ax[0].axhline(y=0, color='r', linestyle='--')
ax[0].set_title("残差 vs 拟合值")
sm.qqplot(residuals, line='s', ax=ax[1])
plt.show()

---

3. 非参数检验（Mann-Whitney U 检验）

场景：比较两组非正态分布数据的差异

from scipy.stats import mannwhitneyu

# 生成非正态数据
group1 = np.random.exponential(2, 50)
group2 = np.random.exponential(3, 50)

# 执行检验
stat, p = mannwhitneyu(group1, group2)
print(f"统计量: {stat}, p 值: {p:.4f}")

---

4. 时间序列分析（ARIMA 模型）

from statsmodels.tsa.arima.model import ARIMA
import pmdarima as pm  # 需要安装 pmdarima

# 生成时间序列数据（模拟股票价格）
np.random.seed(123)
n = 100
dates = pd.date_range(start='2020-01-01', periods=n, freq='D')
data = np.cumsum(np.random.randn(n)) + 50
ts = pd.Series(data, index=dates)

# 自动选择 ARIMA 参数 (p, d, q)
model = pm.auto_arima(ts, seasonal=False, trace=True)

# 拟合模型并预测
model.fit(ts)
forecast = model.predict(n_periods=5)
print("未来5天预测值:\n", forecast)

---

5. 生存分析（Kaplan-Meier 曲线）

from lifelines import KaplanMeierFitter  # 需要安装 lifelines

# 创建模拟生存数据（时间+是否事件发生）
durations = [5, 8, 12, 15, 20, 3, 7, 9]
events = [1, 1, 0, 1, 1, 0, 1, 1]  # 1: 事件发生，0: 删失

# 拟合 Kaplan-Meier 曲线
kmf = KaplanMeierFitter()
kmf.fit(durations, event_observed=events)

# 绘制生存曲线
kmf.plot_survival_function()
plt.title("Kaplan-Meier 生存曲线")
plt.ylabel("生存概率")
plt.show()

---

6. 贝叶斯线性回归（PyMC3 实现）

import pymc3 as pm  # 需要安装 pymc3

# 生成数据
X = np.linspace(0, 10, 50)
y = 3 * X + np.random.normal(0, 2, 50)

# 构建贝叶斯模型
with pm.Model() as model:
    # 定义先验
    alpha = pm.Normal('alpha', mu=0, sigma=10)
    beta = pm.Normal('beta', mu=0, sigma=10)
    sigma = pm.HalfNormal('sigma', sigma=1)
    
    # 定义似然
    mu = alpha + beta * X
    likelihood = pm.Normal('y', mu=mu, sigma=sigma, observed=y)
    
    # 采样
    trace = pm.sample(1000, tune=1000)

# 查看后验分布
pm.plot_posterior(trace, var_names=['alpha', 'beta'])
plt.show()

---

7. 机器学习模型（随机森林回归）

from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error

# 加载波士顿房价数据集
from sklearn.datasets import load_boston
boston = load_boston()
X, y = boston.data, boston.target

# 划分数据集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)

# 训练模型
model = RandomForestRegressor(n_estimators=100)
model.fit(X_train, y_train)

# 预测与评估
preds = model.predict(X_test)
mse = mean_squared_error(y_test, preds)
print(f"测试集 MSE: {mse:.2f}")

---

8. Bootstrap 重抽样

def bootstrap_ci(data, n_bootstrap=1000, ci=95):
    """计算均值的 Bootstrap 置信区间"""
    means = []
    for _ in range(n_bootstrap):
        sample = np.random.choice(data, size=len(data), replace=True)
        means.append(np.mean(sample))
    lower = np.percentile(means, (100 - ci) / 2)
    upper = np.percentile(means, ci + (100 - ci) / 2)
    return lower, upper

# 示例：计算中位数的 95% CI
data = np.random.randn(50) + 2.5
print("Bootstrap 置信区间:", bootstrap_ci(data))

---

运行准备

安装所需库：

pip install numpy scipy pandas statsmodels matplotlib pmdarima lifelines pymc3 scikit-learn

---

总结

通过以上代码，可以快速掌握以下技能：

1. 统计推断：参数估计、假设检验、Bootstrap
2. 建模分析：线性回归、ARIMA、生存分析、贝叶斯模型
3. 机器学习：随机森林集成方法
4. 可视化：残差诊断、生存曲线、后验分布

练习建议：

• 修改代码中的参数（如样本量、分布类型）观察结果变化
• 替换为自己的数据集测试模型
• 查阅官方文档深入理解参数含义（如 statsmodels.OLS、pm.auto_arima）