Model comparison in PyMC3
<p>I am new to PyMC3 and am trying to implement the hierarchical model from Kruschke (2015) section 12.2.2 (model comparison). </p>

<p>I succeeded in defining the full model and then looking at the differences of posterior parameter values (determine whether difference can credibly be said to be zero).</p>

<p>I also tried to explicitly do the comparison in the model as shown in the book (defining a full model and a restricted model and sampling these using a categorical distribution). </p>

<p>Basically I try to implement the below JAGS model definition in PyMC3.<BR>
<a href="http://nbviewer.jupyter.org/github/JWarmenhoven/DBDA-python/blob/master/Notebooks/Chapter%2012.ipynb" rel="nofollow">http://nbviewer.jupyter.org/github/JWarmenhoven/DBDA-python/blob/master/Notebooks/Chapter%2012.ipynb</a><BR>
But I don't know how I can use the model index to select the (pseudo) priors. Any pointers?</p>

<p>JAGS:</p>

<pre><code>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] &lt;-       ( equals(mdlIdx,1)*omega[j] 
                      + equals(mdlIdx,2)*omega0  )   * (kappa[j]-2)+1
    bBeta[j] &lt;- ( 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] &lt;- kappaMinusTwo[j] + 2
    kappaMinusTwo[j] ~ dgamma( 2.618 , 0.0809 ) # mode 20 , sd 20
  }
  # Constants for prior and pseudoprior:
  aP &lt;- 1
  bP &lt;- 1
  # a0[model] and b0[model]
  a0[1] &lt;- .48*500       # pseudo
  b0[1] &lt;- (1-.48)*500   # pseudo 
  a0[2] &lt;- aP            # true
  b0[2] &lt;- bP            # true
  # a[condition,model] and b[condition,model]
  a[1,1] &lt;- aP           # true
  a[2,1] &lt;- aP           # true
  a[3,1] &lt;- aP           # true
  a[4,1] &lt;- aP           # true
  b[1,1] &lt;- bP           # true
  b[2,1] &lt;- bP           # true
  b[3,1] &lt;- bP           # true
  b[4,1] &lt;- bP           # true
  a[1,2] &lt;- .40*125      # pseudo
  a[2,2] &lt;- .50*125      # pseudo
  a[3,2] &lt;- .51*125      # pseudo
  b[1,2] &lt;- (1-.40)*125  # pseudo
  b[2,2] &lt;- (1-.50)*125  # pseudo
  b[3,2] &lt;- (1-.51)*125  # pseudo
  b[4,2] &lt;- (1-.52)*125  # pseudo
  # Prior on model index:
  mdlIdx ~ dcat( modelProb[] )
  modelProb[1] &lt;- .5
  modelProb[2] &lt;- .5
}
</code></pre>

<p>PyMC3:</p>

<pre><code>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.
</code></pre>

<p>Output:</p>

<pre><code>---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
&lt;ipython-input-40-74e77ccc6ce9&gt; in &lt;module&gt;()
      8 
      9     # omega0
---&gt; 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
</code></pre>

<p><strong>UPDATED</strong><br>
After correcting the pseudopriors (missing parenthesis) the results look much better. However, I am not sure whether the pmc.Beta() function works well with arrays as arguments for a and b.
<a href="http://nbviewer.jupyter.org/github/JWarmenhoven/DBDA-python/blob/master/Notebooks/Chapter%2012.ipynb" rel="nofollow">http://nbviewer.jupyter.org/github/JWarmenhoven/DBDA-python/blob/master/Notebooks/Chapter%2012.ipynb</a></p>

<p>The error you are getting is because you are trying to index a list using a tensor. One way to solve this will be to turn the list into a tensor. </p>

<pre><code>import theano.tensor as tt
a0 = tt.as_tensor([.48*500, aP])
</code></pre>

<p>Alternatively you can use <code>pmc.switch()</code> to choose between priors and pseudopriors, something like:</p>

<pre><code>a0 = pm.switch(m_idx, .48*500, aP)
</code></pre>

<p>I did not check your code thoroughly, but notice you have</p>

<pre><code>pmc.switch(eq(m_idx, 0)....)
</code></pre>

<p>Instead, you should write </p>

<pre><code>pmc.switch(pmc.eq(m_idx, 0)....)
</code></pre>

<p>or may be:</p>

<pre><code>pmc.switch(m_idx)....)
</code></pre>

<p>Since 0 evaluates as <code>False</code> and 1 evaluates as <code>True</code>.</p>

<p>Also you have</p>

<pre><code>omega = pmc.Beta('omega0'...)
</code></pre>

<p>And you should have </p>

<pre><code>omega = pmc.Beta('omega'...)
</code></pre>

<p>Your question made me realize I forgot to <a href="https://github.com/aloctavodia/Doing_bayesian_data_analysis" rel="nofollow">port</a> a pseudoprior example. I will do it ASAP.</p>

<p><strong>EDITED</strong></p>

<p>Here if the full model</p>

<pre><code>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)
</code></pre>

<p>Notice that your code had several issues, like missing parenthesis in the definition of the variable <code>b</code> and the order of the prior and pseudopriors was inverted. Additionally I change the code in ordet to let <code>aBeta</code>, <code>bBeta</code> and <code>theta</code>have shape=nCond, and then in the likellihood define <code>p</code> as <code>p=theta[cond_idx]</code>. </p>

<p>I did not check the results against Kruschke's book, but the trace look reasonable.</p>

python bayesian pymc3