我有一些观测数据想要估计参数，我认为这将是尝试 PYMC3 的好机会。
我的数据由一系列记录构成。每条记录包含一对与固定的一小时时间段相关的观察结果。一项观察是给定一小时内发生的事件总数。另一个观察结果是该时间段内的成功数量。例如，一个数据点可能指定在给定的 1 小时内，总共有 1000 个事件，而这 1000 个事件中，有 100 个是成功的。在另一个时间段内，总共可能有1000000个事件，其中120000个是成功的。观察结果的方差不是恒定的，取决于事件的总数，我想控制和建模的部分原因是这种效应。
我这样做的第一步是估计潜在的成功率。我准备了下面的代码，旨在通过使用 scipy 生成两组“观察到的”数据来模拟这种情况。但是，它无法正常工作。
我期望它找到的是：
相反，模型收敛得非常快，但却得出了意想不到的答案。
跟踪图（由于声誉低于 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)
除了任何贝叶斯 MCMC 分析的陷阱之外，您的方法没有任何根本性的错误：(1) 不收敛，(2) 先验，(3) 模型。
不收敛：我找到了一个如下所示的跟踪图：
这不是一件好事，为了更清楚地了解原因，我将更改 Traceplot 代码以仅显示跟踪的后半部分，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 and 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)
我确信还有其他方法在其他研究社区中也会显得更熟悉。
Here is 收集这些评论的笔记本.