我是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
更新
校正伪优先级(缺少括号)后,结果看起来要好得多。但是,我不确定pmc.Beta()函数是否可以将数组用作a和b的参数。http://nbviewer.jupyter.org/github/JWarmenhoven/DBDA-python/blob/master/Notebooks/Chapter%2012.ipynb
您得到的错误是因为您试图使用张量索引列表。解决此问题的一种方法是将列表转换为张量。
import theano.tensor as tt
a0 = tt.as_tensor([.48*500, aP])
或者,您可以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)....)
由于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
并且先验和伪先验的顺序被颠倒了。另外,我将ordet中的代码更改为let aBeta
,bBeta
并theta
具有shape = nCond,然后在可能性中将其定义p
为p=theta[cond_idx]
。
我没有核对克鲁斯克的书的结果,但是踪迹看起来很合理。
本文收集自互联网,转载请注明来源。
如有侵权,请联系[email protected] 删除。
我来说两句