PyMC3中的模型比较

Question

我是PyMC3的新手，正在尝试实现Kruschke（2015）第12.2.2节（模型比较）中的分层模型。在

我成功地定义了完整的模型，然后查看了后验参数值的差异（确定差异是否可以可信地说为零）。在

我还尝试在书中所示的模型中显式地进行比较（定义一个完整的模型和一个受限的模型，并使用分类分布对它们进行抽样）。在

基本上，我尝试在PyMC3中实现以下JAGS模型定义。
http://nbviewer.jupyter.org/github/JWarmenhoven/DBDA-python/blob/master/Notebooks/Chapter%2012.ipynb
但我不知道如何使用模型索引来选择（伪）优先级。有什么建议吗？在

JAGS公司：

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:

^{pr2}$

输出：

---------------------------------------------------------------------------
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

更新
在纠正伪先验（缺少括号）之后，结果看起来好多了。但是，我不确定pmc测试版（）函数可以很好地将数组用作a和b的参数。 http://nbviewer.jupyter.org/github/JWarmenhoven/DBDA-python/blob/master/Notebooks/Chapter%2012.ipynb

Answer 1

您得到的错误是因为您试图使用张量索引列表。解决这个问题的一种方法是将列表转换为张量。在

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

或者，您可以使用pmc.switch()在优先级和伪优先级之间进行选择，类似于：

^{pr2}$

我没有彻底检查你的代码，但是注意你已经检查过了

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

相反，你应该写信

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

或者可能是：

pmc.switch(m_idx)....)

因为0的计算结果为False，而1的计算结果为True。在

你也有

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

你应该有

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

你的问题让我意识到我忘记了一个伪先验的例子。我会尽快做的。在

如果是完整的模型

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)

请注意，您的代码有几个问题，比如变量b的定义中缺少括号，并且prior和pseudo prior的顺序颠倒了。另外，我更改了ordet中的代码，让aBeta、bBeta和{}具有shape=nCond，然后在likelihood中将p定义为p=theta[cond_idx]。在

我没有对照克鲁什克的书检查结果，但痕迹看起来是合理的。在