如何使用PyMC3实现基于威布尔分布回归的简单生存分析

1. 修正模拟数据与公式错误 #

你当前的模拟逻辑不符合生存分析的常规数据结构：生成的y是威布尔概率密度值加噪声，而生存分析的因变量是失效时间，不是密度值。同时你手写的威布尔PDF公式有误，将分子的$(x-\mu)$错写为了$\gamma - \mu$，这也是你后续计算积分出错的核心原因之一。

正确的威布尔失效时间模拟代码如下：

import numpy as np
import pymc3 as pm
import arviz as az

# 威布尔分布真实参数
n = 1000
alpha_true = 1 
gam_true = 0.5
mu_true = 0

# 生成威布尔分布的失效时间
t = mu_true + alpha_true * np.random.weibull(gam_true, size=n)
# 设置右删失规则：失效时间大于7.5的样本做删失处理，仅观测到阈值时间
cens = (t < 7.5).astype(int) # 1代表观测到完整失效时间，0代表右删失
t_obs = np.where(cens == 1, t, 7.5)

2. 威布尔生存函数的解析形式 #

三参数威布尔的CDF和生存函数都有标准解析解，不需要数值积分：

• 累积分布函数：$F(t) = 1 - \exp\left( -\left( \frac{t - \mu}{\alpha} \right)^\gamma \right), \quad t > \mu$
• 生存函数（即你需要的$P(T > t)$，右删失样本的似然计算基础）：$S(t) = 1 - F(t) = \exp\left( -\left( \frac{t - \mu}{\alpha} \right)^\gamma \right)$

3. 带删失的自定义似然实现 #

PyMC3中可以用pm.Potential实现混合似然：非删失样本用概率密度PDF计算对数似然，删失样本用生存函数计算对数似然，实现代码如下：

with pm.Model() as weibull_survival_model:
    # 先验定义
    alpha = pm.Normal('alpha', mu=1, sigma=1)
    gamma = pm.Normal('gamma', mu=0.5, sigma=1)
    mu = pm.Normal('mu', mu=0, sigma=1)
    
    # 拆分非删失与删失样本
    uncens_mask = cens == 1
    t_uncens = t_obs[uncens_mask]
    t_cens = t_obs[~uncens_mask]
    
    # 非删失样本对数似然：log(PDF(t))
    logp_uncens = pm.math.log(gamma/alpha) + (gamma - 1)*pm.math.log((t_uncens - mu)/alpha) - ((t_uncens - mu)/alpha)**gamma
    pm.Potential('uncens_likelihood', pm.math.sum(logp_uncens))
    
    # 删失样本对数似然：log(S(t))
    logp_cens = - ((t_cens - mu)/alpha)**gamma
    pm.Potential('cens_likelihood', pm.math.sum(logp_cens))
    
    # 采样
    trace = pm.sample(2000, return_inferencedata=True)

4. 结果验证 #

采样完成后执行az.summary(trace)即可查看参数后验分布统计值，正常情况下后验均值会和你设定的真实参数alpha=1、gamma=0.5、mu=0高度接近。
如果你确实需要保留原逻辑中加噪声的y作为观测值，只需要将上述代码中的失效时间替换为你的y，对应调整先验范围即可，删失似然的计算逻辑无需修改。

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