6

我正在尝试改编 Cameron Davidson-Pilon 的 《黑客贝叶斯方法》第 1 章“介绍我们的第一个锤子：PyMC”中的文本消息示例来处理多个观察结果。下面的解决方案似乎有效，但我是 pymc 的新手，我不确定这是在 pymc 中处理多个时间序列观察的好方法。任何建议将不胜感激！

回顾一下来自 Bayesian Methods for Hackers 的文本消息示例，观察结果包括 74 天的文本消息计数，如下图所示。

在此处输入图像描述

本书使用一个开关点参数 (tau) 和两个指数参数 (lambda1 和 lambda2) 对这个过程进行建模，它们分别控制 tau 之前和之后的泊松分布消息计数。对于这个例子，pymc 使用以下代码产生近似的解：tau=45、lambda1=18 和 lambda2=23，这与本书的代码几乎相同：

import numpy as np
import pymc

observation = np.loadtxt( './txtdata.csv' ) #data available at the book's GitHub site
n_days      = observation.size    #number of days
alpha       = 1./20  #assume a mean of 20 messages per day
lambda1     = pymc.Exponential("lambda1", alpha)
lambda2     = pymc.Exponential("lambda2", alpha)
tau         = pymc.DiscreteUniform("tau", lower=0, upper=n_days)

@pymc.deterministic
def lambda_(tau=tau, lambda1=lambda1, lambda2=lambda2):
    a       = np.zeros(n_days)
    a[:tau] = lambda1
    a[tau:] = lambda2
    return a
observation_model  = pymc.Poisson("observation", lambda_, value=observation, observed=True)

model   = pymc.Model([observation_model, tau, lambda1, lambda2])
mcmc    = pymc.MCMC(model)
mcmc.sample(40000, 10000)

print()
print( mcmc.trace('tau')[:].mean() )
print( mcmc.trace('lambda1')[:].mean() )
print( mcmc.trace('lambda2')[:].mean() )

我的问题是：这应该如何适应处理多个观察？

我的解决方案出现在下面，并且似乎正在工作，但我想知道是否有更好的方法来模拟 pymc 中的问题。

首先，我使用 tau=45、lambda1=18 和 lambda2=23 生成五个随机观测值，如下所示：

n_observations = 5
n_days  = 74
alpha   = 1./20
lambda1 = pymc.Exponential("lambda1", alpha)
lambda2 = pymc.Exponential("lambda2", alpha)
tau     = pymc.DiscreteUniform("tau", lower=0, upper=n_days)

@pymc.deterministic
def lambda_single(tau=tau, lambda1=lambda1, lambda2=lambda2):
    a       = np.zeros(n_days)
    a[:tau] = lambda1
    a[tau:] = lambda2
    return a

observation_generator  = pymc.Poisson("observation_generator", lambda_single)
tau.set_value(45)
lambda1.set_value(18)
lambda2.set_value(23)
n_observations = 5
observations = np.array(  [observation_generator.random() for i in range(n_observations)]  )

运行上面的代码会产生一个（5 x 74）“观察”数组，代表五个不同人的数据，例如，超过 74 天，如下图所示。

在此处输入图像描述

下一步是我不确定的部分：这五个观察结果应该如何在 pymc 中建模？这是我所拥有的：

@pymc.deterministic
def lambda_multiple(tau=tau, lambda1=lambda1, lambda2=lambda2):
    a = np.zeros( (n_observations, n_days) )
    a[:, :tau] = lambda1
    a[:, tau:] = lambda2
    return a

observation_model  = pymc.Poisson("observations", lambda_multiple, value=observations, observed=True)

运行这个模型似乎会产生 tau、lambda1 和 lambda2 的预期结果，但我想知道这是否是处理多个观察的合适方法？

4

0 回答 0