wehub-skill-sync
659 lines
16 KiB
Markdown
659 lines
16 KiB
Markdown
---
|
||
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 进行贝叶斯推断
|
||
|
||
---
|
||
|
||
## 何时使用贝叶斯方法
|
||
|
||
**使用贝叶斯方法当**：
|
||
- 你有先验信息需要纳入
|
||
- 你想要直接的概率陈述
|
||
- 样本量较小
|
||
- 模型较复杂（层次模型、缺失数据等）
|
||
- 你希望随着数据到来不断更新分析
|
||
|
||
**频率学派方法可能足够当**：
|
||
- 标准分析且样本量大
|
||
- 无先验信息
|
||
- 计算资源有限
|
||
- 审稿人不熟悉贝叶斯方法
|