PyMC:多个时间序列的观察结果(来自“贝叶斯黑客方法”的文本消息示例的改编)

Tod*_*ddP 5 python time-series bayesian pymc probabilistic-programming

我正在尝试改编Cameron Davidson-Pilon的贝叶斯黑客方法,第1章“介绍我们的第一把锤子:PyMC”中的文本消息示例, 以处理多种观察结果。下面的解决方案似乎正在工作,但是我对pymc并不陌生,我不确定这是否是处理pymc中多个时间序列观测值的好方法。任何建议将不胜感激!

为了重述贝叶斯黑客方法的文本消息示例,观察结果包括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() )
Run Code Online (Sandbox Code Playgroud)

我的问题是:应如何调整它以处理多个观察结果?

我的解决方案显示在下面,并且似乎可以正常工作,但是我想知道是否有更好的方法可以在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)]  )
Run Code Online (Sandbox Code Playgroud)

运行上面的代码将产生一个(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)
Run Code Online (Sandbox Code Playgroud)

运行此模型似乎可以产生tau,lambda1和lambda2的预期结果,但是我想知道这是否是处理多个观测值的适当方法?