Биномиальное распределение с параметром loc в pymc3 - PullRequest
Новая
       65

Биномиальное распределение с параметром loc в pymc3

0 голосов
/ 09 июля 2019

Я хотел бы использовать биномиальное распределение, которое сдвигается на параметр loc (как в scipy) в модели pymc3. E.g.:

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

Но Binomial не допускает параметр сдвига.

Я пытался создать его самостоятельно, следуя различным учебным пособиям на сайте 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 ненаблюдаемая и распределена биномиально с (N-z,p), с известным 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, учитывая другие параметры / неизвестные, которые являются просто биномиальным распределением. Но тогда я бы не стал передавать «наблюдаемые» значения, X-z. Другой способ, о котором я подумал, - сначала определить распределения z и X_0, а затем использовать pm.Deterministic для B. Но детерминистский класс не принимает наблюдаемые значения (я полагаю, поскольку он не знает, как оценить вероятность).

1 Ответ

1 голос
/ 11 июля 2019

Вставить скопированный исходный код с помощью pymc3 и добавить параметр 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) # <---
...