10

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

我的数据结构为一系列记录。每条记录都包含一对与固定的一小时时段相关的观察结果。一项观察是在给定小时内发生的事件总数。另一个观察是该时间段内的成功次数。因此，例如，一个数据点可能会指定在给定的 1 小时内，总共有 1000 个事件，在这 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)
4

2 回答 2

33

除了任何贝叶斯 MCMC 分析的缺陷之外，您的方法没有任何根本性错误：（1）不收敛，（2）先验，（3）模型。

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

包含老化的跟踪图

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

删除了老化的跟踪图

先验：收敛的一个主要挑战是您的先验total_lambda_tau，这是贝叶斯建模中的一个典型陷阱。尽管使用 prior 可能看起来没有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)

我敢肯定，在其他研究社区中，还有其他一些看起来更熟悉的方法。

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

于 2014-06-17T19:21:30.543 回答
1

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

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

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

3.) 您为 total_lambda 选择分层先验（即在这些参数上有两个统一先验的正常）确保链需要很长时间才能收敛，除非您明智地选择开始（如第 2 点中所示）。您基本上为 MCMC 链引入了许多不必要的自由度，以使其迷失。鉴于 total_lambda 必须是非自然的，我将为 total_lambda 在合适的范围内选择一个统一先验（例如，从 0 到观察到的最大值）。

4.) 您使用 Metropolis 采样器。20000 个样本可能还不够。尝试 60000 并丢弃前 20000 作为老化。Metropolis Sampler 可能需要一段时间来调整步长，因此很可能它花费了前 20000 个样本来主要拒绝提案和调整。尝试其他采样器，例如 NUTS。

于 2014-06-17T08:07:27.380 回答