pymc3：多个观测值

11 投票

2 回答

5758 浏览

提问于 2025-04-18 09:56

我有一些观察数据，想要估计一些参数，觉得这是试试PYMC3的好机会。

我的数据是以一系列记录的形式存在。每条记录包含一对与固定一小时时间段相关的观察结果。一个观察结果是这个小时内发生的事件总数，另一个观察结果是这个时间段内的成功次数。举个例子，某个数据点可能显示在某个一小时的时间段内，总共有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，我无法发布traceplot，但它看起来并不有趣——快速收敛，且在与输入数据不对应的数字上有明显的峰值。我很好奇我所采取的方法是否存在根本性的问题。以下代码应该如何修改才能得到正确的/预期的结果？

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)

参数估计数据模拟事件计数成功率模型收敛观察数据 pymc3 贝叶斯建模

2 个回答

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

1.) 我觉得成功的计数（叫做错误？）应该遵循一个二项分布（Binomial），也就是总次数为n，成功的概率为p，而不是泊松分布（Poisson）。

2.) 链条是从哪里开始的呢？如果从一个MAP（最大后验估计）或MLE（最大似然估计）的配置开始会比较合理，除非你使用纯Gibbs采样。否则，链条可能需要很长时间来适应，这可能就是这里发生的情况。

3.) 你选择了一个层次先验（hierarchical prior）来处理total_lambda（也就是用两个均匀分布的参数来构建一个正态分布），这会导致链条需要很长时间才能收敛，除非你聪明地选择起始点（就像第二点提到的）。这样做实际上给MCMC链引入了很多不必要的自由度，导致它可能迷失方向。因为total_lambda必须是非负的，我建议为total_lambda选择一个合适范围内的均匀分布先验（比如从0到观察到的最大值）。

4.) 你使用了Metropolis采样器。20000个样本可能不够。试试60000个样本，并把前20000个当作适应期（burn-in）丢掉。Metropolis采样器可能需要一些时间来调整步长，所以它可能在前20000个样本中主要是在拒绝提案和调整参数。你也可以尝试其他采样器，比如NUTS。

回答于 2025-04-18 由 Python大师

分享举报

你的方法本身没有什么根本性的问题，但在进行贝叶斯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]),
                         ])

这样会生成一个看起来可以接受的追踪图：

顺序扫描Metropolis的追踪图

模型：正如@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)

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

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

回答于 2025-04-18 由 Python大师

分享举报

pymc3：多个观测值

2 个回答

撰写回答