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対用量をプロットするが、私は二項分布のたびに描く異なる得るので、これは明らかに正しくありません：私は（バイオアッセイモデルの用語を使用して）これをやってみました。

これは別の質問ですが、どのようにフィットの良さを評価するのですか？私は "始動する" pymc3のチュートリアルでそのようなことを見つけることができなかった。

ありがとうございました！

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を使用することです：

かもしれませ事後予測を使用して

PPC = pm.sample_ppc（トレース、サンプル= 100、モデル= pmmodel）

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

のようなものチェック（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()