pymc3中带有loc参数的二项分布

Question

我想在 pymc3 模型中使用由参数 loc（如 scipy）移动的二项分布。例如：

with pm.Model() as m1:
    prob = pm.Beta('prob',alpha=2,beta=2)
    x = pm.Binomial('x',n=20,p=prob,loc=5)

但是 Binomial 不允许使用 shift 参数。

我尝试按照 pymc3 网站上的各种教程自己构建它，但没有成功（我是使用 pymc3 和 theano 的初学者）。我最近的尝试（可能非常糟糕）

... 
from scipy.stats import binom

class BinoShift(pm.Discrete):
    def __init__(self, n, p, x, *args, **kwargs):
        super(BinoShift, self).__init__(*args, **kwargs)
        self.n = n
        self.p = p
        self.mode = np.round(n*p)
        self.shift = x

    def logp(self, value):
        n = self.n
        p = self.p
        shift = self.shift
        return binom.logpmf(value,n,p,loc=shift)

背景：我对随机变量 X = X_0 + z 有观察，其中 z 是未观察到的潜在变量，X_0 是未观察到的二项分布with (N-z,p) with N 已知。最终目标是在 p 和 z 上获得后验分布。这几乎对应于具有未观察到的集群分配的混合模型问题。 X \sim \sum_z p(z)(z + Bino(p,N-z))。因此，如果我有带移位参数的二项式分布，我设想的 pymc3 模型看起来像

# generate data; kept simple here, but N and z may actually differ across sample 
size = 500
N = 20
p = 0.7
z = 5

X = np.random.binomial(N-z,p,size=size) + z

with pm.Model() as mixture:
    prob = pm.Beta('prob',alpha=2,beta=2)

    weight = pm.Dirichlet('weight',a=np.array([1]*N))
    comp = [pm.Binomial('X_{}'.format(i),n=N-i,p=prob,loc=i) for i in range(N)]
    like = pm.Mixture('like',w=weight,comp_dists=comp,observed=X)

我尝试将此问题构建到 pymc3 模型中的其他方法包括一个分层模型，最后一行引用 X_0 的分布，给定另一个 paramters/unkowns，这只是二项式分布。但是后来我不知道如何传递 "observed" 值，X-z。我想到的另一种方法是先定义 z 和 X_0 的分布，然后使用 pm.Deterministic 来表示 B。但是确定性 class 不接受观察值（我猜是因为它不知道如何评估可能性）。

Answer 1

Copy-pastedpymc3的源代码并添加了loc参数（更改已标记）：

import numpy as np
import theano.tensor as tt

from pymc3.distributions.dist_math import bound, binomln, logpow
from pymc3.math import tround
from pymc3.theanof import floatX, intX
from pymc3.distributions.distribution import Discrete

class BinoShift(Discrete):

    def __init__(self, n, p, loc, *args, **kwargs): # <---
        super().__init__(*args, **kwargs)
        self.n = n = tt.as_tensor_variable(intX(n))
        self.loc = loc = tt.as_tensor_variable(intX(loc)) # <--- 
        self.p = p = tt.as_tensor_variable(floatX(p))
        self.mode = tt.cast(tround(n * p), self.dtype)


    def logp(self, value):
        n = self.n
        p = self.p
        loc = self.loc # <---

        k = value-loc # <---
        return bound(
            binomln(n, k) + logpow(p, k) + logpow(1 - p, n - k),
            0 <= k, k <= n,
            0 <= p, p <= 1) # <---

pymc3中带有loc参数的二项分布

Binomial distribution with loc parameter in pymc3

python

distribution

mixture-model

pymc3