---
name: bayesian-statistical-analysis
description: 贝叶斯统计分析指南——涵盖贝叶斯 vs 频率学派哲学、贝叶斯定理、先验分布、贝叶斯假设检验、贝叶斯估计、常见贝叶斯分析、层次模型、模型比较、模型检查及报告规范
metadata:
  type: reference
---

# 贝叶斯统计分析

本文档为进行和解读贝叶斯统计分析提供指导，贝叶斯统计为频率学派（经典）统计提供了一种替代框架。

## 贝叶斯 vs. 频率学派哲学

### 根本差异

| 方面 | 频率学派 | 贝叶斯学派 |
|--------|-------------|----------|
| **概率解释** | 事件的长期发生频率 | 信念/不确定性的程度 |
| **参数** | 固定但未知 | 具有分布的随机变量 |
| **推断** | 基于抽样分布 | 基于后验分布 |
| **主要输出** | p 值、置信区间 | 后验概率、可信区间 |
| **先验信息** | 不正式纳入 | 通过先验显式纳入 |
| **假设检验** | 拒绝/不拒绝原假设 | 给定数据下假设的概率 |
| **样本量** | 通常需要最小样本量 | 可处理任意样本量 |
| **解释** | 间接（给定 H₀ 下数据的概率） | 直接（给定数据下假设的概率） |

### 关键问题差异

**频率学派**：「如果原假设为真，观察到如此极端或更极端数据的概率是多少？」

**贝叶斯学派**：「给定观测数据，假设为真的概率是多少？」

贝叶斯问题更加直观，直接回答了研究者想要知道的问题。

---

## 贝叶斯定理

**公式**：
```
P(θ|D) = P(D|θ) × P(θ) / P(D)
```

**用文字表述**：
```
后验 = 似然 × 先验 / 证据
```

其中：
- **θ（theta）**：感兴趣的参数（例如均值差、相关系数）
- **D**：观测数据
- **P(θ|D)**: 后验分布（看到数据后对 θ 的信念）
- **P(D|θ)**: 似然（给定 θ 下数据的概率）
- **P(θ)**: 先验分布（看到数据前对 θ 的信念）
- **P(D)**: 边际似然/证据（归一化常数）

---

## 先验分布

### 先验的类型

#### 1. 信息性先验

**何时使用**：当你拥有大量先验知识时，来源包括：
- 先前的研究
- 专家知识
- 理论
- 预试验数据

**示例**：元分析显示效应量 d ≈ 0.40，SD = 0.15
- 先验：Normal(0.40, 0.15)

**优点**：
- 纳入已有知识
- 更高效（所需样本量更小）
- 可在小样本情况下稳定估计

**缺点**：
- 具有主观性（但主观性也可成为优势）
- 必须被合理证明并保持透明
- 若强先验与数据冲突，可能引发争议

---

#### 2. 弱信息性先验

**何时使用**：大多数应用场景的默认选择

**特征**：
- 正则化估计（防止极端值）
- 在中等样本量下对后验影响极小
- 防止计算问题

**示例先验**：
- 效应量：Normal(0, 1) 或 Cauchy(0, 0.707)
- 方差：Half-Cauchy(0, 1)
- 相关系数：Uniform(-1, 1) 或 Beta(2, 2)

**优点**：
- 在客观性与正则化之间取得平衡
- 计算稳定
- 广泛可接受

---

#### 3. 无信息先验（平坦/均匀先验）

**何时使用**：试图保持「客观」时

**示例**：Uniform(-∞, ∞) 对任意值

**⚠️ 注意**：
- 可能导致非正常后验
- 可能产生不合理的结果
- 并非真正的「无信息」（仍然做了假设）
- 在现代贝叶斯实践中通常不推荐

**更好的替代方案**：使用弱信息性先验

---

### 先验敏感性分析

**始终进行**：检验结果如何随不同先验变化

**流程**：
1. 用默认/计划先验拟合模型
2. 用更分散的先验拟合模型
3. 用更集中的先验拟合模型
4. 比较后验分布

**报告**：
- 若结果相似：证据稳健
- 若结果差异显著：数据不足以压倒先验

**Python 示例**：
```python
import pymc as pm

prior_specs = [
    ('weakly_informative', 0, 1),
    ('diffuse', 0, 10),
    ('informative', 0.5, 0.3),
]

results = {}
for name, mu_prior, sigma_prior in prior_specs:
    with pm.Model() as model:
        effect = pm.Normal('effect', mu=mu_prior, sigma=sigma_prior)
        # ... 似然函数和观测数据
        trace = pm.sample(2000, tune=1000, return_inferencedata=True)
        results[name] = trace
```

---

## 贝叶斯假设检验

### 贝叶斯因子（BF）

**定义**：两个竞争假设的证据比率

**公式**：
```
BF₁₀ = P(D|H₁) / P(D|H₀)
```

**解释**：

| BF₁₀ | 证据强度 |
|------|----------|
| >100 | 支持 H₁ 的决定性证据 |
| 30-100 | 支持 H₁ 的极强证据 |
| 10-30 | 支持 H₁ 的强证据 |
| 3-10 | 支持 H₁ 的中等证据 |
| 1-3 | 支持 H₁ 的微弱证据 |
| 1 | 无证据 |
| 1/3-1 | 支持 H₀ 的微弱证据 |
| 1/10-1/3 | 支持 H₀ 的中等证据 |
| 1/30-1/10 | 支持 H₀ 的强证据 |
| 1/100-1/30 | 支持 H₀ 的极强证据 |
| <1/100 | 支持 H₀ 的决定性证据 |

**相对于 p 值的优势**：
1. 可提供支持原假设的证据
2. 不依赖于抽样意图（无「偷看」问题）
3. 直接量化证据
4. 可随更多数据更新

**Python 计算**：
```python
# Pingouin 0.5+：用于独立双边 t 检验的 BF10；单边 BF 已被移除。
import pingouin as pg

result = pg.ttest(group1, group2, correction=False)
bf10 = result['BF10'].values[0]

# 严格贝叶斯因子：BayesFactor（R）、JASP 或 PyMC 模型比较（参见 pymc 技能）
```

---

### 实际等价区间（ROPE）

**目的**：定义可忽略效应量的范围

**流程**：
1. 定义 ROPE（例如 d ∈ [-0.1, 0.1] 为可忽略效应）
2. 计算后验落在 ROPE 内的百分比
3. 做出判断：
   - >95% 在 ROPE 内：接受实际等价
   - >95% 在 ROPE 外：拒绝等价
   - 其他情况：无定论

**优势**：直接检验实际显著性

**Python 示例**：
```python
# 定义 ROPE
rope_lower, rope_upper = -0.1, 0.1

# 计算后验在 ROPE 内的百分比
in_rope = np.mean((posterior_samples > rope_lower) &
                  (posterior_samples < rope_upper))

print(f"{in_rope*100:.1f}% 的后验在 ROPE 内")
```

---

## 贝叶斯估计

### 可信区间

**定义**：以 X% 的概率包含参数的区间

**95% 可信区间的解释**：
> 「真实参数有 95% 的概率落在此区间内。」

**这正是人们以为置信区间所表示的含义**（但在频率学派框架中并非如此）

**类型**：

#### 等尾区间（ETI）
- 第 2.5 至第 97.5 百分位
- 计算简单
- 对于偏态分布可能不包含众数

#### 最高密度区间（HDI）
- 包含分布 95% 的最窄区间
- 始终包含众数
- 更适合偏态分布

**Python 计算**：
```python
import arviz as az

# 等尾区间
eti = np.percentile(posterior_samples, [2.5, 97.5])

# HDI
hdi = az.hdi(posterior_samples, hdi_prob=0.95)
```

---

### 后验分布

**解读后验分布**：

1. **集中趋势**：
   - 均值：后验的平均值
   - 中位数：第 50 百分位
   - 众数：最可能的值（MAP——最大后验估计）

2. **不确定性**：
   - 标准差：后验的离散程度
   - 可信区间：量化不确定性

3. **形状**：
   - 对称：近似正态
   - 偏态：非对称不确定性
   - 多峰：存在多个合理值

**可视化**：
```python
import matplotlib.pyplot as plt
import arviz as az

# 带 HDI 的后验图
az.plot_posterior(trace, hdi_prob=0.95)

# 迹线图（检查收敛性）
az.plot_trace(trace)

# 森林图（多个参数）
az.plot_forest(trace)
```

---

## 常见贝叶斯分析

### 贝叶斯 T 检验

**目的**：比较两组（t 检验的贝叶斯替代方法）

**输出**：
1. 均值差的后验分布
2. 95% 可信区间
3. 贝叶斯因子（BF₁₀）
4. 方向性假设的概率（例如 P(μ₁ > μ₂)）

**Python 实现**：
```python
import pymc as pm
import arviz as az

# 贝叶斯独立样本 t 检验
with pm.Model() as model:
    # 组均值的先验
    mu1 = pm.Normal('mu1', mu=0, sigma=10)
    mu2 = pm.Normal('mu2', mu=0, sigma=10)

    # 合并标准差的前沿
    sigma = pm.HalfNormal('sigma', sigma=10)

    # 似然函数
    y1 = pm.Normal('y1', mu=mu1, sigma=sigma, observed=group1)
    y2 = pm.Normal('y2', mu=mu2, sigma=sigma, observed=group2)

    # 衍生量：均值差
    diff = pm.Deterministic('diff', mu1 - mu2)

    # 后验抽样
    trace = pm.sample(2000, tune=1000, return_inferencedata=True)

# 分析结果
print(az.summary(trace, var_names=['mu1', 'mu2', 'diff']))

# group1 > group2 的概率
prob_greater = np.mean(trace.posterior['diff'].values > 0)
print(f"P(μ₁ > μ₂) = {prob_greater:.3f}")

# 绘制后验图
az.plot_posterior(trace, var_names=['diff'], ref_val=0)
```

---

### 贝叶斯方差分析（ANOVA）

**目的**：比较三个或更多组

**模型**：
```python
import pymc as pm

with pm.Model() as anova_model:
    # 超先验
    mu_global = pm.Normal('mu_global', mu=0, sigma=10)
    sigma_between = pm.HalfNormal('sigma_between', sigma=5)
    sigma_within = pm.HalfNormal('sigma_within', sigma=5)

    # 组均值（层次结构）
    group_means = pm.Normal('group_means',
                            mu=mu_global,
                            sigma=sigma_between,
                            shape=n_groups)

    # 似然函数
    y = pm.Normal('y',
                  mu=group_means[group_idx],
                  sigma=sigma_within,
                  observed=data)

    trace = pm.sample(2000, tune=1000, return_inferencedata=True)

# 后验对比
contrast_1_2 = trace.posterior['group_means'][:,:,0] - trace.posterior['group_means'][:,:,1]
```

---

### 贝叶斯相关分析

**目的**：估计两个变量之间的相关系数

**优势**：提供相关系数值的分布

**Python 实现**：
```python
import pymc as pm

with pm.Model() as corr_model:
    # 相关系数的先验
    rho = pm.Uniform('rho', lower=-1, upper=1)

    # 转换为协方差矩阵
    cov_matrix = pm.math.stack([[1, rho],
                                [rho, 1]])

    # 似然函数（二元正态分布）
    obs = pm.MvNormal('obs',
                     mu=[0, 0],
                     cov=cov_matrix,
                     observed=np.column_stack([x, y]))

    trace = pm.sample(2000, tune=1000, return_inferencedata=True)

# 相关系数汇总
print(az.summary(trace, var_names=['rho']))

# 相关系数为正的概率
prob_positive = np.mean(trace.posterior['rho'].values > 0)
```

---

### 贝叶斯线性回归

**目的**：建模预测变量与结果变量之间的关系

**优势**：
- 所有参数均带有不确定性
- 自然正则化（通过先验）
- 可纳入先验知识
- 预测的可信区间

**Python 实现**：
```python
import pymc as pm

with pm.Model() as regression_model:
    # 系数的先验
    alpha = pm.Normal('alpha', mu=0, sigma=10)  # 截距
    beta = pm.Normal('beta', mu=0, sigma=10, shape=n_predictors)
    sigma = pm.HalfNormal('sigma', sigma=10)

    # 期望值
    mu = alpha + pm.math.dot(X, beta)

    # 似然函数
    y_obs = pm.Normal('y_obs', mu=mu, sigma=sigma, observed=y)

    trace = pm.sample(2000, tune=1000, return_inferencedata=True)

# 后验预测检验
with regression_model:
    ppc = pm.sample_posterior_predictive(trace)

az.plot_ppc(ppc)

# 带不确定性的预测
with regression_model:
    pm.set_data({'X': X_new})
    posterior_pred = pm.sample_posterior_predictive(trace)
```

---

## 层次（多水平）模型

**何时使用**：
- 嵌套/聚类数据（学生嵌套于学校）
- 重复测量
- 元分析
- 跨组的变异性效应

**关键概念**：部分池化
- 完全池化：忽略组别（有偏）
- 无池化：分别分析各组（高方差）
- 部分池化：跨组借用信息（贝叶斯）

**示例：变截距模型**：
```python
with pm.Model() as hierarchical_model:
    # 超先验
    mu_global = pm.Normal('mu_global', mu=0, sigma=10)
    sigma_between = pm.HalfNormal('sigma_between', sigma=5)
    sigma_within = pm.HalfNormal('sigma_within', sigma=5)

    # 组级截距
    alpha = pm.Normal('alpha',
                     mu=mu_global,
                     sigma=sigma_between,
                     shape=n_groups)

    # 似然函数
    y_obs = pm.Normal('y_obs',
                     mu=alpha[group_idx],
                     sigma=sigma_within,
                     observed=y)

    trace = pm.sample()
```

---

## 模型比较

### 方法

#### 1. 贝叶斯因子
- 直接比较模型证据
- 对先验设定敏感
- 计算量可能较大

#### 2. 信息准则

**WAIC（广泛适用信息准则）**：
- AIC 的贝叶斯类比
- 越小越好
- 考虑了参数的有效数量

**LOO（留一法交叉验证）**：
- 估计样本外预测误差
- 越小越好
- 比 WAIC 更稳健

**Python 计算**：
```python
import arviz as az

# 计算 WAIC 和 LOO
waic = az.waic(trace)
loo = az.loo(trace)

print(f"WAIC: {waic.elpd_waic:.2f}")
print(f"LOO: {loo.elpd_loo:.2f}")

# 比较多个模型
comparison = az.compare({
    'model1': trace1,
    'model2': trace2,
    'model3': trace3
})
print(comparison)
```

---

## 检查贝叶斯模型

### 1. 收敛诊断

**R-hat（Gelman-Rubin 统计量）**：
- 比较链内方差与链间方差
- 接近 1.0 的值表明收敛
- R-hat < 1.01：良好
- R-hat > 1.05：收敛不佳

**有效样本量（ESS）**：
- 独立样本的数量
- 越高越好
- 建议每链 ESS > 400

**迹线图**：
- 应呈现「毛虫状」
- 无趋势，无卡顿链

**Python 检查**：
```python
# 带诊断的自动汇总
print(az.summary(trace, var_names=['parameter']))

# 可视化诊断
az.plot_trace(trace)
az.plot_rank(trace)  # 秩图
```

---

### 2. 后验预测检验

**目的**：模型生成的数据是否与观测数据相似？

**流程**：
1. 从后验生成预测
2. 与实际数据比较
3. 寻找系统性偏差

**Python 实现**：
```python
with model:
    ppc = pm.sample_posterior_predictive(trace)

# 可视化检验
az.plot_ppc(ppc, num_pp_samples=100)

# 定量检验
obs_mean = np.mean(observed_data)
pred_means = [np.mean(sample) for sample in ppc.posterior_predictive['y_obs']]
p_value = np.mean(pred_means >= obs_mean)  # 贝叶斯 p 值
```

---

## 报告贝叶斯结果

### T 检验报告示例

> 「采用贝叶斯独立样本 t 检验比较 A 组和 B 组。使用了弱信息性先验：均值差采用 Normal(0, 1)，合并标准差采用 Half-Cauchy(0, 1)。均值差的后验分布均值为 5.2（95% CI [2.3, 8.1]），表明 A 组得分高于 B 组。贝叶斯因子 BF₁₀ = 23.5 为组间差异提供了强证据，A 组均值超过 B 组均值的概率为 99.7%。」

### 回归报告示例

> 「采用弱信息性先验（系数使用 Normal(0, 10)，残差标准差使用 Half-Cauchy(0, 5)）拟合了贝叶斯线性回归模型。模型解释了可观的方差（R² = 0.47，95% CI [0.38, 0.55]）。学习时长（β = 0.52，95% CI [0.38, 0.66]）和先前的 GPA（β = 0.31，95% CI [0.17, 0.45]）是可信的预测变量（95% CI 排除零）。后验预测检验显示模型拟合良好。收敛诊断令人满意（所有 R-hat < 1.01，ESS > 1000）。」

---

## 优势与局限

### 优势

1. **直观的解释**：关于参数的直接概率陈述
2. **纳入先验知识**：利用所有可用信息
3. **灵活**：轻松处理复杂模型
4. **无 p 值操纵**：可随时查看数据
5. **量化不确定性**：完整的后验分布
6. **小样本**：可处理任意样本量

### 局限

1. **计算量大**：需要 MCMC 抽样（可能较慢）
2. **先验设定**：需要思考和论证
3. **复杂性**：学习曲线较陡
4. **软件工具**：相比频率学派方法工具较少
5. **沟通成本**：可能需要向审稿人/读者做解释

---

## 关键 Python 包

使用 uv 安装（参见 SKILL.md）。ArviZ 0.23+ 需要 Python 3.12+。

- **PyMC** (`pymc>=5`)：完整的贝叶斯建模框架
- **ArviZ** (`arviz>=0.17`)：可视化与诊断工具（[文档](https://python.arviz.org)）
- **Bambi**：回归模型的高级接口（`uv pip install bambi`）
- **PyStan**：Stan 的 Python 接口
- **TensorFlow Probability**：使用 TensorFlow 进行贝叶斯推断

---

## 何时使用贝叶斯方法

**使用贝叶斯方法当**：
- 你有先验信息需要纳入
- 你想要直接的概率陈述
- 样本量较小
- 模型较复杂（层次模型、缺失数据等）
- 你希望随着数据到来不断更新分析

**频率学派方法可能足够当**：
- 标准分析且样本量大
- 无先验信息
- 计算资源有限
- 审稿人不熟悉贝叶斯方法