pymc3：多个观测值

Question

我有一些观测数据，我想估计参数，我认为这将是一个很好的机会来尝试PYMC3。在

我的数据是由一系列记录构成的。每个记录包含一对与固定的一小时周期相关的观察结果。一个观察是给定时间内发生的事件总数。另一个观察是在这段时间内成功的次数。例如，一个数据点可以指定在给定的1小时内，总共有1000个事件，而在这1000个事件中，100个是成功的。在另一个时间段内，总共可能有1000000个事件，其中120000个是成功的。观测值的方差不是恒定的，它取决于事件的总数，我想控制和建模的部分原因就是这种影响。在

我这样做的第一步是估计潜在的成功率。我准备了下面的代码，通过使用scipy生成两组“观察到的”数据来模拟这种情况。但是，它不能正常工作。
我希望它找到的是：

• 损失λ系数约为0.1
• 总λ约为120。在

相反，模型收敛得很快，但得到了一个意想不到的答案。在

• 总λ和总λμ分别是5e5附近的尖峰。在
• 损失λ系数约为0。在

traceplot（由于声誉低于10，我无法发布）相当无趣-快速收敛，在与输入数据不符的数字处出现尖峰。我很好奇我采取的方法是否根本上是错误的。应如何修改以下代码以获得正确/预期的结果？在

from pymc import Model, Uniform, Normal, Poisson, Metropolis, traceplot 
from pymc import sample 
import scipy.stats

totalRates = scipy.stats.norm(loc=120, scale=20).rvs(size=10000)
totalCounts = scipy.stats.poisson.rvs(mu=totalRates) 
successRate = 0.1*totalRates 
successCounts = scipy.stats.poisson.rvs(mu=successRate) 

with Model() as success_model: 
    total_lambda_tau= Uniform('total_lambda_tau', lower=0, upper=100000)
    total_lambda_mu = Uniform('total_lambda_mu', lower=0, upper=1000000)
    total_lambda = Normal('total_lambda', mu=total_lambda_mu, tau=total_lambda_tau)
    total = Poisson('total', mu=total_lambda, observed=totalCounts) 

    loss_lambda_factor = Uniform('loss_lambda_factor', lower=0, upper=1)
    success_rate = Poisson('success_rate', mu=total_lambda*loss_lambda_factor, observed=successCounts) 

with success_model: 
    step =  Metropolis() 
    success_samples = sample(20000, step) #, start)


plt.figure(figsize=(10, 10)) 
_ = traceplot(success_samples)

Answer 1

你的方法没有本质上的错误，除了任何贝叶斯MCMC分析的缺陷：（1）不收敛，（2）先验，（3）模型。在

不收敛：我发现了一个如下所示的跟踪打印：

这不是一件好事，为了更清楚地了解原因，我将更改traceplot代码以仅显示跟踪的后半部分，traceplot(success_samples[10000:])：

先验知识：收敛的一个主要挑战是你在total_lambda_tau上的先验知识，这是贝叶斯建模中的一个典型陷阱。虽然使用previor Uniform('total_lambda_tau', lower=0, upper=100000)可能显得不具信息性，但其效果是说您非常确定total_lambda_tau很大。例如，它小于10的概率是0.0001。更改之前的

total_lambda_tau= Uniform('total_lambda_tau', lower=0, upper=100)
total_lambda_mu = Uniform('total_lambda_mu', lower=0, upper=1000)

结果得到一个更有希望的追踪图：

然而，这仍然不是我在traceplot中寻找的，为了得到更满意的结果，我建议使用“顺序扫描大都市”步骤（这是PyMC2在类似模型中默认的步骤）。可以按如下方式指定：

^{pr2}$

这将生成一个似乎可以接受的跟踪打印：

模型：正如@KaiLondenberg所回应的那样，您在total_lambda_tau和{}上采用的方法不是标准方法。您描述的事件总数变化很大（一小时1000个，下一个小时1000000个），但是您的模型假定它是正态分布的。在空间流行病学中，我所看到的类似数据的方法是一个更像这样的模型：

import pymc as pm, theano.tensor as T
with Model() as success_model: 
    loss_lambda_rate = pm.Flat('loss_lambda_rate')
    error = Poisson('error', mu=totalCounts*T.exp(loss_lambda_rate), 
            observed=successCounts)

我相信，在其他研究领域，还有其他方法似乎更为熟悉。在

这是a notebook collecting up these comments。在

Answer 2

我发现这个模型有几个潜在的问题。在

1）我认为成功是重要的（称为错误？）应遵循二项式分布（n=总计，p=损失λ系数），而不是泊松分布。在

2.） 链条从哪里开始？除非使用纯Gibbs采样，否则从MAP或MLE配置开始是有意义的。否则，链条可能需要很长时间老化，这可能就是这里发生的情况。在

3.）您对total_lambda的层次优先权的选择（即在这些参数上有两个统一的优先权）可以确保链将花费很长的时间来收敛，除非您明智地选择开始（如第2点）。你实际上引入了很多不必要的自由度，让MCMC链迷失方向。考虑到total_lambda必须是非ngative的，我将在适当的范围内（例如从0到观察到的最大值）为total_lambda选择一个统一的先验。在

4.）你用大都会采样器。2万个样本可能不够。尝试60000，然后丢弃前20000个作为烧坏。Metropolis Sampler可能需要一段时间来调整步长，因此很可能它花费了前20000个样本来主要拒绝建议和调整。试试其他的采样器，比如坚果。在