如何修复PyMC中截断正态分布的模型参数化？贝叶斯元分析场景

问题概述 #

我正在构建用于事件报告率的贝叶斯元分析模型，核心目标是估计每个研究中参与者报告事件X的概率。每个参与者需提供20-40个二元响应（X或非X），因此每个研究可计算出受试者的X报告平均率，我需要搭建模型来估计该平均率。

当前建模逻辑 #

• 概率分布需约束在0-1区间，匹配报告率的取值范围。
• 由于报告率通常集中在0-0.1区间，使用对数转换后的形式比原始概率更适合建模。
• 采用delta方法近似对数标准误，公式为SE/p。

现有模拟数据与模型代码 #

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

import pymc as pm
import arviz as az

np.random.seed(123)
study_means = np.exp(np.random.normal(np.log(.1),.1,20))
study_se = np.random.uniform(0.02, 0.05, 20)

def leave_one_out(df):
    
    for i in range(len(df)):
        # 返回删除第i行后的新DataFrame
        yield df.drop(index=df.index[i]).copy()


def fit_meta_model(df):

    coords = {"random_effects": df.index.values}

    with pm.Model(coords=coords) as reparam_model:
    
        # 数据定义
        log_rate = pm.Data("log_rate", df["log_rate"].values, dims="random_effects")
        log_se = pm.Data("log_se", df["log_se"].values, dims="random_effects")
        
        # 先验分布
        theta_p = pm.Beta("theta_p", alpha=1, beta=1) 
        theta = pm.Deterministic("theta", pm.math.log(theta_p))
                       
        tau = pm.HalfCauchy("tau", beta=.5)
    
        # 非中心化参数化
        z = pm.Normal("z", mu=0, sigma=1, dims="random_effects")
        theta_k_raw = theta + z * tau
        theta_k = pm.Deterministic("theta_k", pm.math.minimum(theta_k_raw, 0.0), dims="random_effects")
        
        # 似然函数（截断正态分布）
        obs = pm.TruncatedNormal("obs", mu=theta_k, sigma=log_se, upper=0.0, observed=log_rate, dims="random_effects")
        
        # 采样
        prior = pm.sample_prior_predictive(500)
        idata = pm.sample(1000, tune=1000, target_accept=.95)

    return idata, prior

# 创建模拟数据集并加入极端值
study_means = np.exp(np.random.normal(np.log(.1),.1,18))
study_means = np.r_[study_means,.6,.6]

study_se = np.random.uniform(0.02, 0.1,size=20)

df = pd.DataFrame({
    "mean_rate":study_means,
    "se":study_se,
})

df["log_rate"] = np.log(df["mean_rate"])
df["log_se"] = df["se"] / df["mean_rate"]

# 拟合模型
idata,prior = fit_meta_model(df) 

# 绘制先验分布
f,a = plt.subplots(1,3,figsize=(15,5))

az.plot_posterior(np.exp(prior.prior.theta),ax=a[0])
az.plot_posterior(prior.prior.tau,ax=a[1])

az.plot_posterior(np.exp(prior.prior_predictive.obs),ax=a[2])

当前模型痛点 #

• 潜在参数theta通过Beta(1,1)（均匀分布）先验转换为对数形式，tau采用标准HalfCauchy(0, .5)建模，两者先验符合预期。
• 非中心化参数化和似然的截断正态分布实现存在难度，部分模拟数据集会出现采样发散，提高target_accept可缓解该问题。
• 先验预测模拟结果存在偏斜，希望改用负对数正态分布（更贴合报告率集中在0-0.1的分布特征）但尚未实现，询问是否可通过重新参数化或模型修改来降低先验预测的偏斜程度。

解决方案 #

1. 直接采用负对数正态分布建模原始率

既然报告率集中在0-0.1区间，负对数正态分布是更合适的选择，可以直接对原始率建模，避免对数转换带来的偏斜问题：

def fit_meta_model_reparam(df):
    coords = {"random_effects": df.index.values}
    
    with pm.Model(coords=coords) as model:
        # 数据
        mean_rate = pm.Data("mean_rate", df["mean_rate"].values, dims="random_effects")
        se = pm.Data("se", df["se"].values, dims="random_effects")
        
        # 群体水平先验：负对数正态的参数
        mu_log = pm.Normal("mu_log", mu=np.log(0.1), sigma=0.5)
        sigma_log = pm.HalfCauchy("sigma_log", beta=0.5)
        
        # 研究水平参数：每个研究的对数率
        log_theta_k = pm.Normal("log_theta_k", mu=mu_log, sigma=sigma_log, dims="random_effects")
        theta_k = pm.Deterministic("theta_k", pm.math.exp(log_theta_k), dims="random_effects")
        
        # 似然函数：原始率的正态近似（样本量足够大时中心极限定理适用）
        obs = pm.Normal("obs", mu=theta_k, sigma=se, observed=mean_rate, dims="random_effects")
        
        prior = pm.sample_prior_predictive(500)
        idata = pm.sample(1000, tune=1000, target_accept=.95)
    
    return idata, prior

该模型直接对原始率的对数建模，利用负对数正态的特性（取值范围0到正无穷，且集中在低值区），先验预测会更贴合实际数据分布。

2. 优化原对数域模型的先验和截断逻辑

如果坚持在对数域建模，可调整先验并优化截断方式：

• 将theta_p的先验改为Beta(2, 20)，让原始率的先验更集中在0-0.1区间，避免均匀分布带来的极端值影响。
• 替换pm.math.minimum为更平滑的截断，直接用pm.TruncatedNormal生成theta_k：

theta_k = pm.TruncatedNormal("theta_k", mu=theta, sigma=tau, upper=0.0, dims="random_effects")

这样能避免生硬截断带来的采样问题和先验偏斜。

3. 先验预测检查的调整

在对数域做先验预测时，需将结果转换回原始率再评估，同时可以对先验参数设置更贴合业务认知的约束，比如给tau设置上限，避免研究间变异过大导致的偏斜。

内容的提问来源于stack exchange，提问作者Paradeisios