了解 PyMC 中的因子潜力

机器算法验证 时间序列 贝叶斯 泊松分布 变化点 pymc
2022-04-12 14:58:05

我正在尝试从 PyMC 文档中了解因子潜力，但在实现部分需要一些帮助 - 否则可能会证明我误解了潜力是如何工作的。

想象一下，我们正在构建PyMC 文档教程中指定的泊松开关点模型。现在我们要引入一个与文档中描述的完全一样的因子势，使得早期和晚期泊松均值之间的差异小于 1。

当我尝试在下面的代码中实现潜力时，我看到允许约束之外的值出现在后验分布中。为什么 potentialCheck 的后验分布包含这些值？显然，我做错了什么...... 

from pymc import DiscreteUniform, Exponential, deterministic, Poisson, Uniform, potential
import numpy as np

disasters_array =   \
 np.array([ 4, 5, 4, 0, 1, 4, 3, 4, 0, 6, 3, 3, 4, 0, 2, 6,
               3, 3, 5, 4, 5, 3, 1, 4, 4, 1, 5, 5, 3, 4, 2, 5,
               2, 2, 3, 4, 2, 1, 3, 2, 2, 1, 1, 1, 1, 3, 0, 0,
               1, 0, 1, 1, 0, 0, 3, 1, 0, 3, 2, 2, 0, 1, 1, 1,
               0, 1, 0, 1, 0, 0, 0, 2, 1, 0, 0, 0, 1, 1, 0, 2,
               3, 3, 1, 1, 2, 1, 1, 1, 1, 2, 4, 2, 0, 0, 1, 4,
               0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1])

switchpoint = DiscreteUniform('switchpoint', lower=0, upper=110, doc='Switchpoint[year]')

early_mean = Exponential('early_mean', beta=1.)
late_mean = Exponential('late_mean', beta=1.)

@deterministic(plot=False)
def rate(s=switchpoint, e=early_mean, l=late_mean):
    ''' Concatenate Poisson means '''
    out = empty(len(disasters_array))
    out[:s] = e
    out[s:] = l
    return out

@potential(plot=True)
def examplePotential(em=early_mean, lm=late_mean ):
    if abs(em-lm) < 1:
         return e
    else:
         return 1.

@deterministic(plot=True)
def potentialCheck(e=early_mean, l=late_mean):
    '''Replicate the Potential to plot and check if the constraint has held.'''
    return abs(e-l)

disasters = Poisson('disasters', mu=rate, value=disasters_array, observed=True)

disaster_model = [switchpoint, early_mean, late_mean, examplePotential, potentialCheck]

from pymc import MCMC
M = MCMC(disaster_model)
M.sample(iter=10000, burn=1000, thin=10)

from pymc.Matplot import plot as mcplot

mcplot(M)
3个回答

这些天我都在学习相同的教程，所以我远非专家，但我在您的代码中看到了一个可能的错字：

@potential(plot=True)
def examplePotential(em=early_mean, lm=late_mean ):
    if abs(em-lm) < 1:
         return e  <<<<<<< return em !!!!!!!! 
    else:
         return 1.

有了这个改变，它对我来说很好。

因子潜力应返回对数概率，而不是参数值。请记住，它只是一个任意的对数概率项——它本身没有价值。

例如，如果您想将特定值限制为大于零，您可以有一个因子潜力，它为参数的所有非正值值返回 -inf。

我同意克里斯的回答。如果您对代码进行以下更改，它应该可以工作：

@potential(plot=True)
def examplePotential(em=early_mean, lm=late_mean ):
    if abs(em-lm) < 1:
        return 0. ####### No penalty applied if conditions are met
    else:
        return -np.inf ##### Infinite potential applied if not met

您还必须更改 early_mean 和 late_mean 的初始值，以使它们满足潜在条件。我只是将这两个数量的值都设为 1。

其它你可能感兴趣的问题
上一篇为什么正态分布变量的抽样分布自动也是正态分布 下一篇了解 ACF - ARIMA 模型的 PCF 图