I have a list of n observations, each of which is the sum of two Weibull-distributed variables:

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设置或两者都有问题 . 我尝试使用电位调整this advice来模拟潜在变量 . 首先,我为每个观察定义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)

You can see my full code here. 我通过模拟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不包括真值 . 正如this histogram所示,由此产生的分布甚至没有在正确的范围内 . (蓝色显示我使用的模拟数据x [i],而红色显示MCMC运行中代表性迭代的完全不同的推断分布 . )

使用不同的随机种子再次运行,我得到shape1 = 4.65,scale1 = 23.3,shape2 = 0.83,scale2 = 21.3 . 这种分布有点接近事实 . 有没有办法改变MCMC设置,以便始终为这类问题获得不错的结果?任何关于更有效地使用PyMC的建议都非常感谢 .

更新 - 尝试了“辅助”MCMC运行:

我还尝试通过使用接近事实的值初始化人口级参数来协助MCMC运行 . 结果稍好一些,但我现在发现了系统的偏见 . 下面的直方图显示了观测值(蓝色)与拟合分布(红色)的真实分布 . 右尾很合适,但是合身无法捕捉左侧的尖峰 . 对于n = 60和100的种群大小,这种偏差一致地发生 . 我不确定这是更多的PyMC问题还是一般的MCMC算法问题 .

Histogram of "assisted" run