pymc3:多个观察值

Question

我有一些观察数据,我想估计参数,我认为这是一个尝试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.

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代码以仅显示跟踪的后半部分traceplot(success_samples[10000:]):

先前:融合的一个主要挑战是你的先验total_lambda_tau,这是贝叶斯建模的一个典范陷阱.虽然使用先前可能看起来很无法提供信息Uniform('total_lambda_tau', lower=0, upper=100000),但这样做的效果就是说你很确定它total_lambda_tau很大.例如,它小于10的概率是.0001.改变之前的

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

导致更有希望的traceplot:

然而,这仍然不是我在traceplot中寻找的东西,并且为了获得更令人满意的东西,我建议使用"顺序扫描Metropolis"步骤(这是PyMC2默认用于类似模型的步骤).您可以按如下方式指定:

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

这会产生一个似乎可以接受的traceplot:

该模型:作为@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)

我相信还有其他方法在其他研究社区中也会更为熟悉.

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