Я новичок в PyMC3 и пытаюсь реализовать иерархическую модель от Kruschke (2015), раздел 12.2.2 (сравнение модели).Сравнение моделей в PyMC3
Мне удалось определить полную модель, а затем посмотреть на различия значений задних параметров (определить, можно ли достоверно сказать, может ли разница быть нулевой).
Я также попытался явно провести сравнение в модели, как показано в книге (определение полной модели и ограниченной модели и выборка их с использованием категориального распределения).
В основном я пытаюсь реализовать определение модели JAGS ниже в PyMC3.
http://nbviewer.jupyter.org/github/JWarmenhoven/DBDA-python/blob/master/Notebooks/Chapter%2012.ipynb
Но я не знаю, как я могу использовать индекс модели, чтобы выбрать (pseudo) priors. Любые указатели?
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:
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 (отсутствует скобка) результаты выглядят гораздо Бетт р. Однако я не уверен, работает ли функция pmc.Beta() с массивами как аргументы для a и b. http://nbviewer.jupyter.org/github/JWarmenhoven/DBDA-python/blob/master/Notebooks/Chapter%2012.ipynb
Спасибо, определяя эти приоритеты как тензорные константы, которые скомпилировали модель. (Я довольно новичок в деталях Теано). Но при индексировании этих констант все еще что-то не так, используя модельный индекс 'm_idx'. Эта переменная не изменяет значение при выборке. – Jordi
Я определил приоритеты, как вы предположили: вложение предварительных констант a, b, a0, b0 с помощью 'tt.as_tensor()'. Остальная модель не изменилась. – Jordi
Yikes ... Я полностью пропускаю эти ошибки! Будет проверяться как можно скорее и отчитываться. – Jordi