Как выполнить зачетные результаты для дискретных значений в pymc3?

Question

Я совершенно не знаком с pymc3, поэтому, пожалуйста, извините, что это, вероятно, тривиально. У меня очень простая модель, где я предсказываю функцию двоичного ответа. Модель представляет собой почти дословную копию этого примера: https://github.com/pymc-devs/pymc3/blob/master/pymc3/examples/gelman_bioassay.py

Я возвращаю параметры модели (альфа, бета и тета), но я не могу понять, как переопределить предсказания модели против. входные данные. Я попытался сделать это (используя жаргон модели биотестирования):

from scipy.stats import binom 

mean_alpha = mean(trace['alpha']) 
mean_beta = mean(trace['beta']) 

pred_death = binom.rvs(n, 1./(1.+np.exp(-(mean_alpha + mean_beta * dose)))) 

, а затем черчения дозы против pred_death, но это явно не правильно, как я получаю разные розыгрыши биномиального распределения каждый раз.

Связанный с этим другой вопрос, как я могу оценить доброту соответствия? Кажется, я не мог найти ничего подобного в учебнике pimc3, начатом с самого начала.

Большое спасибо за любой совет!

Answer 1

Если вы хотите построить дозу против pred_death, где pred_death вычисляется из средних оценочных значений альфа и бета, а затем сделать:

pred_death = 1./(1. + np.exp(-(mean_alpha + mean_beta * dose))) 
plt.plot(dose, pred_death) 

вместо этого, если вы хотите построить дозу против pred_death, где pred_death является рассчитанный с учетом неопределенности в заднем для альфа и бета. Тогда, вероятно, самый простой способ заключается в использовании функции sample_ppc:

Может быть что-то вроде

ррс = pm.sample_ppc (след, образцы = 100, модель = pmmodel)

for i in range(100): 
    plt.plot(dose, ppc['deaths'][i], 'bo', alpha=0.5) 

Использование POSTERIOR Predictive Проверки (ppc) - это способ проверить, насколько хорошо ваша модель ведет себя, сравнивая предсказания модели с вашими фактическими данными. Here у вас есть пример sample_ppc

Другие варианты могут заключаться в построении среднего значения плюс некоторый интервал интереса.

Answer 2

Привет простой способ сделать это состоит в следующем:

from pymc3 import * 
from numpy import ones, array 

# Samples for each dose level 
n = 5 * ones(4, dtype=int) 
# Log-dose 
dose = array([-.86, -.3, -.05, .73]) 


def invlogit(x): 
    return np.exp(x)/(1 + np.exp(x)) 



with Model() as model: 

    # Logit-linear model parameters 
    alpha = Normal('alpha', 0, 0.01) 
    beta = Normal('beta', 0, 0.01) 

    # Calculate probabilities of death 
    theta = Deterministic('theta', invlogit(alpha + beta * dose)) 

    # Data likelihood 
    deaths = Binomial('deaths', n=n, p=theta, observed=[0, 1, 3, 5]) 
    start = find_MAP() 
    step = NUTS(scaling=start) 
    trace = sample(2000, step, start=start, progressbar=True) 




import matplotlib.pyplot as plt 

death_fit = np.percentile(trace.theta,50,axis=0) 

plt.plot(dose, death_fit,'g', marker='.', lw='1.25', ls='-', ms=5, mew=1) 

plt.show()