Как определить пользовательский приоритет в PyMC3

Я хотел бы знать, можно ли определить пользовательское значение в PyMC3 (и как это сделать). Из здесь кажется, что в PyMC2 это относительно легко сделать (без необходимости изменять исходный код code), но в PyMC3 не все так просто (или я чего-то не понимаю). Я пытаюсь воспроизвести априор из книги «Выполнение байесовского анализа данных», которая реализована в ОШИБКАХ:

model {
# Likelihood. Each flip is Bernoulli.
for ( i in 1 : N1 ) { y1[i]  ̃ dbern( theta1 ) }
for ( i in 1 : N2 ) { y2[i]  ̃ dbern( theta2 ) }
# Prior. Curved scallo not ps!
x  ̃ dunif(0,1)
y  ̃ dunif(0,1)
N <- 4
xt <- sin( 2*3.141593*N * x ) / (2*3.141593*N) + x
yt <- 3 * y + (1/3)
xtt <- pow( xt , yt )
theta1 <- xtt
theta2 <- y
}

Априор не имеет особого смысла, это просто пример того, как определить собственный априор и универсальность ОШИБОК.

Моя попытка реализовать вышеуказанный пользовательский пример:

from __future__ import division
import numpy as np
import pymc as pm
from pymc import Continuous
from theano.tensor import sin, log

# Generate the data
y1 = np.array([1, 1, 1, 1, 1, 0, 0])  # 5 heads and 2 tails
y2 = np.array([1, 1, 0, 0, 0, 0, 0])  # 2 heads and 5 tails

class Custom_prior(Continuous): 
"""
custom prior
"""
    def __init__(self, y, *args, **kwargs):
        super(Custom_prior, self).__init__(*args, **kwargs)
        self.y = y
        self.N = 4
        self.mean = 0.625  # FIXME
    def logp(self, value):
        N = self.N
        y = self.y
        return -log((sin(2*3.141593*N * value)
                     / (2*3.141593*N) + value)**(3 * y + (1/3)))

with pm.Model() as model:
    theta2 = pm.Uniform('theta2', 0, 1)  # prior
    theta1 = Custom_prior('theta1', theta2)  # prior
    # define the likelihood
    y1 = pm.Bernoulli('y1', p=theta1, observed=y1)
    y2 = pm.Bernoulli('y2', p=theta2, observed=y2)
    # Generate a MCMC chain
    start = pm.find_MAP()  # Find starting value by optimization
    trace = pm.sample(5000, pm.NUTS(), progressbar=False)

РЕДАКТИРОВАТЬ

После ответа chris-fonnesbeck

Я думаю, мне нужно что-то вроде:

with pm.Model() as model:
    theta2 = pm.Uniform('theta2', 0, 1)  # prior
    N = 4
    theta1 = pm.DensityDist('theta1', lambda value: -log((sin(2*3.141593*N * value)
                       / (2*3.141593*N) + value)**(3 * theta2 + (1/3))))
    # define the likelihood
    y1 = pm.Bernoulli('y1', p=theta1, observed=y1)
    y2 = pm.Bernoulli('y2', p=theta2, observed=y2)

    # Generate a MCMC chain
    start = pm.find_MAP()  # Find starting value by optimization
    trace = pm.sample(10000, pm.NUTS(), progressbar=False) # Use NUTS sampling

Единственная проблема заключается в том, что я получаю одно и то же значение для всех апостериорных выборок тета1 и тета2, я предполагаю, что есть какая-то проблема с моим пользовательским априором или комбинацией априорных значений и вероятностей. Успешное определение пользовательского априора можно найти в этом пример


python bayesian pymc3
person aloctavodia    schedule 12.07.2014    source источник

Ответы (1)


arrow_upward
4
arrow_downward

Можете ли вы опубликовать полную модель BUGS? В приведенном выше это просто выглядит как серия детерминированных преобразований в BUGS после априорных значений для x и y, а не определение априорного значения.

Предполагая, что logp выше — это то, что вам нужно, вы можете реализовать его в PyMC гораздо проще:

def logp(value, y):
    N  = 4
    return -log((sin(2*3.141593*N * value)
                 / (2*3.141593*N) + value)**(3 * y + (1/3)))

theta1 = pm.DensityDist('theta1', logp, value, y=theta2)
person Chris Fonnesbeck    schedule 14.07.2014
comment
Большое спасибо! Я редактирую свой вопрос, чтобы добавить полную модель. В любом случае, как указано в книге, априорное значение не является реальным априорным значением, которое можно было бы использовать в исследованиях, это просто пример того, что априорное распределение может быть любым, а не обязательно хорошо известным распределением, таким как бета, гауссово и т. д. - person aloctavodia; 14.07.2014