pymc3: 多个观测值

Question

我有一些观测数据需要估计参数，想着试试PYMC3。我的数据结构是一系列记录。每个记录包含一对与固定一小时期间相关的观测值。一个观测值是在给定小时内发生的事件总数，另一个观测值是该时间段内成功的数量。例如，一个数据点可能指定在给定的一个1小时周期内，总共有1000个事件，其中100个是成功的。在另一个时间段，总共可能有1000000个事件，其中120000个是成功的。观测值的方差不是恒定的，并且取决于总事件数，部分原因是我想要控制和建模的效应。

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

• loss_lambda_factor大致为0.1
• total_lambda（和total_lambda_mu）大致为120。

相反，模型非常快地收敛，但却得到了意外的答案：

• total_lambda和total_lambda_mu分别是5e5左右的尖峰。
• loss_lambda_factor大约为0。

跟踪图（由于声望低于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)

- sozen

Answer 1

除了任何贝叶斯MCMC分析的缺陷之外，您的方法基本上没有什么问题：（1）非收敛，（2）先验，（3）模型。

非收敛：我发现跟踪图看起来像这样：

这不是一件好事，为了更清楚地看到原因，我会改变追踪图代码，只显示追踪的后半部分，traceplot(success_samples[10000:])：

先验分布: 收敛的一个主要挑战是您对total_lambda_tau的先验分布，这是贝叶斯建模中的一个典型陷阱。尽管使用先验分布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)

结果是一个更有前途的追踪图：

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

step =  pm.CompoundStep([pm.Metropolis([total_lambda_mu]),
                         pm.Metropolis([total_lambda_tau]),
                         pm.Metropolis([total_lambda]),
                         pm.Metropolis([loss_lambda_factor]),
                         ]) 

这将生成一个看起来合理的跟踪图：

这个模型: 正如@KaiLondenberg所回答的那样，你对total_lambda_tau和total_lambda_mu使用先验概率的方法并不是标准的。你描述了事件总数差异很大（一个小时为1,000次，下一个为1,000,000次），但你的模型假设其服从正态分布。在空间流行病学中，我看到的类似数据的处理方式更像是这样的模型:

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)

我相信在其他研究领域也有更为熟悉的方法。

这里是一个收集了这些评论的笔记本。

- Abraham D Flaxman

---

亚伯拉罕，你的回答非常详细，太感谢了！ - sozen

Answer 2

我看到这个模型存在几个潜在问题。

1.) 我认为成功次数（称为误差？）应该遵循二项式分布Binomial(n=total,p=loss_lambda_factor)，而不是泊松分布。

2.) 链的起点在哪里？从MAP或MLE配置开始是有意义的，除非您使用纯Gibbs抽样。否则，链可能需要很长时间进行燃烧，这可能就是这里发生的情况。

3.) 您选择使用层次先验(total_lambda)（即具有两个参数均匀先验的正态分布），这确保了链需要很长时间才能收敛，除非您明智地选择起始点（如第2点）。您实际上引入了许多不必要的自由度，使MCMC链容易迷失方向。考虑到total_lambda必须是非负的，我建议在适当的范围内（例如从0到观察到的最大值）选择Uniform先验。

4.) 你使用了大都市采样器。20000个样本可能不足以满足需求。尝试使用60000个样本，并将前20000个作为烧掉的样本进行丢弃。大都市采样器可能需要一段时间来调整步长，因此它很可能会花费前20000个样本主要用于拒绝提议和调整。尝试使用其他采样器，如NUTS。

- Kai Londenberg

---

非常感谢你，凯！ - sozen

"尝试其他采样器，比如NUTS" - 我认为NUTS或HMC不适用于离散分布（泊松分布或二项式分布）。 - merv