0

下面是定义两个随机伯努利随机变量的一种方法，一个依赖于另一个带有装饰器的变量。该模型旨在：

p(A) = 0.5
p(B=True|A=True) = 0.75
p(B=True|A=False) = 0.05

在 PyMC 中使用装饰器是：

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

def make_model():
    @pymc.stochastic(dtype=bool)
    def A(value=False):
        def logp(value):
            return pymc.bernoulli_like(value, 0.5)

        def random():
            return bool(np.random.binomial(1, 10**logp(value)))


    @pymc.stochastic(dtype=bool)
    def B(value=False, A=False):
        def logp(value, A):
            if A:
                logp = pymc.bernoulli_like(value, 0.75)
            else:
                logp = pymc.bernoulli_like(value, 0.05)
            return logp

        def random(A):
            return bool(np.random.binomial(1, 10**logp(True, A)))

    return locals()

test = pymc.Model(make_model())

这是最正确和最紧凑的方法吗？是否可以通过将节点显式定义为pymc.Bernoulli变量而不是使用stochastic装饰器来保存一些代码？

上面的代码似乎是多余的，必须始终定义如何从伯努利 RV 中采样，但这也许是不可避免的？

最后是否有必要从numpylike调用采样函数，np.random.binomial或者 PyMC 中是否有采样函数？

4

1 回答 1

1

我不认为这是正确的。例如，如果我使用以下代码从您的模型中采样，我发现样本的平均值A不是 0.5：

In [20]: pymc.MCMC(test).sample(10000)
 [-----------------100%-----------------] 10000 of 10000 complete in 0.6 sec
In [22]: test.A.trace().mean()
Out[22]:
0.19670000000000001

这是实现模型的一种紧凑方法：

import pymc as pm

def make_model():
    A = pm.Bernoulli('A', .5)
    B_given_A_true = pm.Bernoulli('B_given_A_true', .75)
    B_given_A_false = pm.Bernoulli('B_given_A_false', .05)

    @deterministic
    def B(A=A, B_given_A_true=B_given_A_true,
          B_given_A_false=B_given_A_false):
        if A:
            return B_given_A_true
        else:
            return B_given_A_false

    return locals()

test = pymc.Model(make_model())
于 2015-03-03T04:55:31.453 回答