用可变速率拟合上限泊松过程答案

【问题标题】：Fitting a capped Poisson process with a variable rate用可变速率拟合上限泊松过程
【发布时间】：2014-05-05 20:41:13
【问题描述】：

我正在尝试使用最大后验估计来估计泊松过程的速率，其中速率随时间变化。这是一个速率线性变化的简化示例（λ = ax+b）：

import numpy as np
import pymc

# Observation
a_actual = 1.3
b_actual = 2.0
t = np.arange(10)
obs = np.random.poisson(a_actual * t + b_actual)

# Model
a = pymc.Uniform(name='a', value=1., lower=0, upper=10)
b = pymc.Uniform(name='b', value=1., lower=0, upper=10)


@pymc.deterministic
def linear(a=a, b=b):
    return a * t + b

r = pymc.Poisson(mu=linear, name='r', value=obs, observed=True)

model = pymc.Model([a, b, r])
map = pymc.MAP(model)
map.fit()
map.revert_to_max()

print "a :", a._value
print "b :", b._value

这工作正常。但是我的实际泊松过程受到确定性值的限制。由于我无法将观察到的值与确定性函数相关联，因此我正在为我的观察添加一个方差很小的正态随机函数：

import numpy as np
import pymc

# Observation
a_actual = 1.3
b_actual = 2.0
t = np.arange(10)
obs = np.random.poisson(a_actual * t + b_actual).clip(0, 10)

# Model
a = pymc.Uniform(name='a', value=1., lower=0, upper=10)
b = pymc.Uniform(name='b', value=1., lower=0, upper=10)


@pymc.deterministic
def linear(a=a, b=b):
    return a * t + b

r = pymc.Poisson(mu=linear, name='r')


@pymc.deterministic
def clip(r=r):
    return r.clip(0, 10)

rc = pymc.Normal(mu=r, tau=0.001, name='rc', value=obs, observed=True)

model = pymc.Model([a, b, r, rc])
map = pymc.MAP(model)
map.fit()
map.revert_to_max()

print "a :", a._value
print "b :", b._value

此代码产生以下错误：

Traceback (most recent call last):
  File "pymc-bug-2.py", line 59, in <module>
    map.revert_to_max()
  File "pymc/NormalApproximation.py", line 486, in revert_to_max
    self._set_stochastics([self.mu[s] for s in self.stochastics])
  File "pymc/NormalApproximation.py", line 58, in __getitem__
    tot_len += self.owner.stochastic_len[p]
KeyError: 0

知道我做错了什么吗？

【问题讨论】：

标签： pymc

【解决方案1】：

“封顶”是指截断泊松吗？看来这就是你所说的。如果它是左截断（更常见），您可以使用TruncatedPoisson 分布，但由于您正在执行右截断，所以不能（我们应该让它更通用！）。您正在尝试的将不起作用 - Poisson 对象没有 clip() 方法。您可以做的是使用因子潜力。它看起来像这样：

@pymc.potential
def clip(r=r):
    if np.any(r>10):
        return -np.inf
    return 0

这会将r 的值限制为小于10。有关Potential 类的信息，请参阅pymc docs。

【讨论】：