PyMC:估计总体参数,其中每个观测值是两个Weibull分布变量的总和

Question

我有一个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,并且我运行1e5次AdaptiveMetropolis迭代.结果收敛于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算法问题.