5

我有一个 n 观察的列表，每个观察都是两个 Weibull 分布变量的总和：

x[i] = t1[i] + t2[i]
t1[i] ~ Weibull(shape1, scale1)
t2[i] ~ Weibull(shape2, scale2)

我的目标是：

1）估计两个威布尔分布的形状和尺度参数（shape1，scale1，shape2，scale2），

2) 对于每个观测值 x[i]，估计 t1[i]（并且由此得出 t2[i]）。

（旁白：每次观察x[i]是癌症诊断的年龄，t1[i]和t2[i]是肿瘤发展的两个不同时间段。实际模型也涉及到突变数据，但在我之前试试看，我想确保我可以使用 PyMC 来解决这个更简单的问题。）

我正在使用 PyMC2 进行这些估计，看起来运行收敛，但结果不正确。我不知道我的 PyMC 模型语法、MCMC 设置或两者都有问题。我尝试将这个建议改编为使用潜力对潜在变量进行建模。首先，我为每个观察定义 x[i] 和 t1[i]：

for i in xrange(n):
    x[i] = pm.Index('x_%i'%i, x=data, index=i) # data is a list of observations
    t1[i] = pm.Weibull('t1_%i'%i, alpha=shape1, beta=scale1)
    # Ensure that initial guess for t1 is not more than the observed sum:
    if t1[i].value >= x[i].value:
        t1[i].value = 0.95 * x[i].value

然后我为 t2[i] = x[i] - t1[i] 定义一个确定性：

for i in xrange(n):
    def subtractfunc(t1=t1, x=x, ii=i):
        return x[ii] - t1[ii]
    t2[i] = pm.Lambda('t2_%i'%i, subtractfunc)

最后我定义了 t2[i] 的潜力：

t2dist = np.empty(n, dtype=object)
for i in xrange(n):
    def weibfunc(t2=t2, shape2=shape2, scale2=scale2, ii=i):
        return pm.weibull_like(t2[ii], alpha=shape2, beta=scale2)
    t2dist[i] = pm.Potential(logp = weibfunc,
                               name = 't2dist_%i'%i,
                               parents = {'shape2':shape2, 'scale2':scale2, 't2':t2},
                               doc = 'weibull potential for t2',
                               verbose = 0,
                               cache_depth = 2)

你可以在这里看到我的完整代码。我通过模拟 60 个独立观测值进行测试，shape1 = 1，scale1 = 30，shape2 = 6.5，scale2 = 10，并且我运行 AdaptiveMetropolis 的 1e5 次迭代。结果收敛到 shape1=1.94、scale1=37.9、shape2=0.55、scale2=36.1 的平均值，并且 95% HPD 不包括真实值。如该直方图所示，此结果分布甚至不在正确的范围内。（蓝色显示我使用的模拟数据 x[i]，而红色显示与 MCMC 运行中的代表性迭代完全不同的推断分布。）

使用不同的随机种子再次运行，我得到 shape1=4.65，scale1=23.3，shape2=0.83，scale2=21.3。这种分布在某种程度上更接近事实。是否有某种方法可以更改 MCMC 设置以始终为此类问题获得不错的结果？非常感谢任何有关更有效地使用 PyMC 的建议。

更新——尝试了“辅助”MCMC 运行：

我还尝试通过使用接近真实值的值初始化人口级参数来帮助 MCMC 运行。结果稍微好一些，但我现在发现了系统偏差。下面的直方图显示了观测值（蓝色）与拟合分布（红色）的真实分布。右尾非常适合，但无法捕捉到左侧的尖峰。对于人口规模 n = 60 和 100，这种偏差始终存在。我不确定这更多是 PyMC 问题还是一般 MCMC 算法问题。

直方图

4

0 回答 0