PyMC3でのモデル比較

Question

私はPyMC3を初めて使用しており、Kruschke（2015）のセクション12.2.2（モデル比較）から階層モデルを実装しようとしています。

私は完全なモデルを定義し、事後パラメータ値の差異（差異が確実にゼロであると言えるかどうかを判断する）を調べることに成功しました。

私はまた、本のモデル（完全なモデルと制限されたモデルを定義し、これらをカテゴリ型の分布を使ってサンプリングする）のようにモデル内で明示的に比較しようとしました。

基本的に私はPyMC3で以下のJAGSモデル定義を実装しようとしています。
http://nbviewer.jupyter.org/github/JWarmenhoven/DBDA-python/blob/master/Notebooks/Chapter%2012.ipynb
しかし、モデルインデックスを使用して（擬似）プリオを選択する方法はわかりません。すべてのポインタ？

ぎざぎざ：

model { 
    for (s in 1:nSubj) { 
    nCorrOfSubj[s] ~ dbin(theta[s] , nTrlOfSubj[s]) 
    theta[s] ~ dbeta(aBeta[CondOfSubj[s]] , bBeta[CondOfSubj[s]]) 
    } 

    for (j in 1:nCond) { 
    # Use omega[j] for model index 1, omega0 for model index 2: 
    aBeta[j] <-  (equals(mdlIdx,1)*omega[j] 
         + equals(mdlIdx,2)*omega0 ) * (kappa[j]-2)+1 
    bBeta[j] <- (1 - (equals(mdlIdx,1)*omega[j] 
         + equals(mdlIdx,2)*omega0 )) * (kappa[j]-2)+1 
    omega[j] ~ dbeta(a[j,mdlIdx] , b[j,mdlIdx]) 
    } 

    omega0 ~ dbeta(a0[mdlIdx] , b0[mdlIdx]) 
    for (j in 1:nCond) { 
    kappa[j] <- kappaMinusTwo[j] + 2 
    kappaMinusTwo[j] ~ dgamma(2.618 , 0.0809) # mode 20 , sd 20 
    } 
    # Constants for prior and pseudoprior: 
    aP <- 1 
    bP <- 1 
    # a0[model] and b0[model] 
    a0[1] <- .48*500  # pseudo 
    b0[1] <- (1-.48)*500 # pseudo 
    a0[2] <- aP   # true 
    b0[2] <- bP   # true 
    # a[condition,model] and b[condition,model] 
    a[1,1] <- aP   # true 
    a[2,1] <- aP   # true 
    a[3,1] <- aP   # true 
    a[4,1] <- aP   # true 
    b[1,1] <- bP   # true 
    b[2,1] <- bP   # true 
    b[3,1] <- bP   # true 
    b[4,1] <- bP   # true 
    a[1,2] <- .40*125  # pseudo 
    a[2,2] <- .50*125  # pseudo 
    a[3,2] <- .51*125  # pseudo 
    b[1,2] <- (1-.40)*125 # pseudo 
    b[2,2] <- (1-.50)*125 # pseudo 
    b[3,2] <- (1-.51)*125 # pseudo 
    b[4,2] <- (1-.52)*125 # pseudo 
    # Prior on model index: 
    mdlIdx ~ dcat(modelProb[]) 
    modelProb[1] <- .5 
    modelProb[2] <- .5 
} 

がPyMC3：

with pmc.Model() as model_1: 
    # constants 
    aP, bP = 1, 1 

    # Pseudo- and true hyperpriors per model 
    a0 = [.48*500, aP]  
    b0 = [(1-.48)*500, bP] 

    # Lower level pseudo- and true priors per model/condition combination 
    a = np.c_[np.tile(aP, 4), [(.40*125), (.50*125), (.51*125), (.52*125)]] 
    b = np.c_[np.tile(bP, 4), [(1-.40)*125, (1-.50)*125, (1-.51)*125, (1-.52)*125]] 

    # Prior on model index [0,1] 
    m_idx = pmc.Categorical('m_idx', np.asarray([.5, .5])) 

    # Priors on concentration parameters 
    kappa = pmc.Gamma('kappa', 2.618, 0.0809, shape=nCond) 

    # omega0 
    omega0 = pmc.Beta('omega0', a0[m_idx], b0[m_idx])  

    # omega (condition specific) 
    omega = pmc.Beta('omega', a[:,m_idx], b[:,m_idx], shape=nCond) 

    # theta 
    aBeta = pmc.switch(eq(m_idx, 0), omega0 * kappa[cond_idx]+1, omega[cond_idx] * kappa[cond_idx]+1) 
    bBeta = pmc.switch(eq(m_idx, 0), (1-omega0) * kappa[cond_idx]+1, (1-omega[cond_idx]) * kappa[cond_idx]+1) 

    theta = pmc.Beta('theta', aBeta[cond_idx], bBeta[cond_idx], shape=df.index.size) 

    # Likelihood 
    y = pmc.Binomial('y', n=df.nTrlOfSubj.values, p=theta, observed=df.nCorrOfSubj)  




Applied log-transform to kappa and added transformed kappa_log_ to model. 

出力：

--------------------------------------------------------------------------- 
TypeError         Traceback (most recent call last) 
<ipython-input-40-74e77ccc6ce9> in <module>() 
     8 
     9  # omega0 
---> 10  omega0 = pmc.Beta('omega0', a0[m_idx], b0[m_idx]) 
    11 
    12  # omega (condition specific) 

TypeError: list indices must be integers or slices, not FreeRV 

（括弧の欠落）pseudopriorsを修正した後
を更新された結果がはるかにベティを見てr。しかし、pmc.Beta（）関数がaとbの引数として配列でうまく動作するかどうかはわかりません。 http://nbviewer.jupyter.org/github/JWarmenhoven/DBDA-python/blob/master/Notebooks/Chapter%2012.ipynb

Answer 1

エラーは、テンソルを使用してリストのインデックスを作成しようとしているためです。これを解決する1つの方法は、リストをテンソルに変えることです。

import theano.tensor as tt 
a0 = tt.as_tensor([.48*500, aP]) 

また、あなたが何かのように、事前確率とpseudopriorsの間で選択するpmc.switch()を使用することができます：私は徹底的にあなたのコードをチェックしませんでした

a0 = pm.switch(m_idx, .48*500, aP) 

ていますが、その代わり

pmc.switch(eq(m_idx, 0)....) 

を持って気づきます、書く必要があります

pmc.switch(pmc.eq(m_idx, 0)....) 

であってもよい。

pmc.switch(m_idx)....) 

Falseと1と0が評価さTrueとして評価するため。

はまた、あなたは

omega = pmc.Beta('omega0'...) 

を持っていて、

omega = pmc.Beta('omega'...) 

あなたの質問は私がport pseudoprior例に忘れて実感しましたが必要です。私はできるだけ早くそれをやります。ここで

EDITED

あなたのコードは、変数bの定義に括弧が欠落などのいくつかの問題を、持っていたし、前とpseudopriorsの順であったことがフルモデル

with pmc.Model() as model_1: 

# constants 
aP, bP = 1., 1. 

# Pseudo- and true hyperpriors per model 
a0 = tt.as_tensor([aP, .48*500])  
b0 = tt.as_tensor([bP, (1-.48)*500]) 

# Lower level pseudo- and true priors per model/condition combination 
a = tt.as_tensor(np.c_[[(.40*125), (.50*125), (.51*125), (.52*125)], np.tile(aP, 4)]) 
b = tt.as_tensor(np.c_[[((1-.40)*125), ((1-.50)*125), ((1-.51)*125), ((1-.52)*125)], np.tile(bP, 4)]) 

# Prior on model index [0,1] 
m_idx = pmc.Categorical('m_idx', p=np.array([.5, .5])) 

# Priors on concentration parameters 
kappa = pmc.Gamma('kappa', 2.618, 0.0809, shape=nCond) 

# omega0 
omega0 = pmc.Beta('omega0', a0[m_idx], b0[m_idx])  

# omega (condition specific) 
omega = pmc.Beta('omega', a[:,m_idx], b[:,m_idx], shape=nCond) 

# theta 
aBeta = pmc.switch(pmc.eq(m_idx, 0), omega0 * kappa+1, omega * kappa+1) 
bBeta = pmc.switch(pmc.eq(m_idx, 0), (1-omega0) * kappa+1, (1-omega) * kappa+1) 

theta = pmc.Beta('theta', aBeta, bBeta, shape=nCond) 

# Likelihood 
y = pmc.Binomial('y', n=df.nTrlOfSubj.values, p=theta[cond_idx], observed=df.nCorrOfSubj) 
trace = pmc.sample(1000) 

気付いた場合反転した。さらに、aBeta,bBetaおよびthetaにshape = nCondを設定し、pをp=theta[cond_idx]と定義するように、ordetのコードを変更します。

私はKruschkeの本との結果をチェックしませんでしたが、トレースは妥当なものです。