Проверка модели PyMC Bernoulli

Question

В настоящее время я пытаюсь выполнить проверку модели с помощью PyMC, где моя модель - модель Бернулли, и у меня есть бета-версия. Я хочу сделать и (i) график gof, а также (ii) рассчитать заднее предсказательное значение p.

У меня есть мой код, работающий с моделью Binomial, но я изо всех сил стараюсь найти правильный способ создания модели Bernoulli. К сожалению, нет ни одного примера, где бы я ни работал. Мой код выглядит следующим образом:

import pymc as mc 
import numpy as np 
alpha = 2 
beta = 2 
n = 13 
yes = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,0,0]) 

p = mc.Beta('p',alpha,beta) 
surv = mc.Bernoulli('surv',p=p,observed=True,value=yes) 
surv_sim = mc.Bernoulli('surv_sim',p=p) 

mc_est = mc.MCMC({'p':p,'surv':surv,'surv_sim':surv_sim}) 
mc_est.sample(10000,5000,2) 

import matplotlib.pylab as plt 
plt.hist(mc_est.surv_sim.trace(),bins=range(0,3),normed=True) 
plt.figure() 
plt.hist(mc_est.p.trace(),bins=100,normed=True) 

mc.Matplot.gof_plot(mc_est.surv_sim.trace(), 10/13., name='surv') 

#here I have issues 
D = mc.discrepancy(yes, surv_sim, p.trace()) 
mc.Matplot.discrepancy_plot(D) 

Основная проблема, которую я имею в определении ожидаемых значений функции discrepancy. Просто использование p.trace() здесь не работает, так как это вероятности. Так или иначе, мне нужно включить размер выборки здесь, но я изо всех сил стараюсь сделать это так же, как и для модели Binomial. Я также не совсем уверен, если я правильно делаю gof_plot.

Надеюсь, кто-то может помочь мне здесь! Благодаря!

Answer 1

согласно discrepancy функции doc string, параметры:

observed : Iterable of observed values (size=(n,)) 
simulated : Iterable of simulated values (size=(r,n)) 
expected : Iterable of expected values (size=(r,) or (r,n)) 

Таким образом, вы должны исправить две вещи:

1) изменить ваши результаты моделирования, чтобы иметь размер n (например, 13 в вашем примере):

surv_sim = mc.Bernoulli('surv_sim', p=p, size=n)

2) инкапсулироватьс bernoulli_expval методом:

D = mc.discrepancy(yes, surv_sim.trace(), mc.bernoulli_expval(p.trace()))

(. bernoulli_expval просто плюет назад p)

С этими двумя пересадками, я получаю следующее: