Покрытие кода OpenBUGS в PYMC3

Question

Ранее я работал с OpenBUGS/WinBUGS для выполнения моей байесовской статистики, но решил перейти к использованию пакета PYMC3 в Python. Так что я довольно новичок в pacakage и все еще учусь использовать его полностью. У меня возникли некоторые трудности с преобразованием моего кода ошибок в PYMC3. Оригинальный код ошибки:

model {
for (i in 1 : N) {
    Tobs[i] ~ dpois(mu[i])
    mu[i]<- u[i]*lamda[i]
    u[i]~dbern(p[i])
    log(lamda[i]) <- a0 + a1*insta[i] + a2*shear[i] + b[i]
    p[i]<-exp(-beta/exp(popd[i]))}
b[1:N] ~ car.normal(adj[], weights[], num[], tau)
for(k in 1:sumNumNeigh) {weights[k] <- 1}
a0 ~ dnorm(0, 0.001)
a1 ~ dnorm(0, 0.001)
a2 ~ dnorm(0, 0.001)
beta<-exp(betaaux)
betaaux~dnorm(1, 0.001)
tau ~ dgamma(0.01,0.01)
sigma <-sqrt(1 / tau)
}

Я написал это в Python как так:

model1 = pm.Model()
with model1:
    #Vague prior
        betaaux = pm.Normal('betaaux', mu=1.0, tau=1.0e-3)
        beta = pm.Deterministic('beta', tt.exp(betaaux))
    #Random effects (hierarchial) prior
        tau_c = pm.InverseGamma('tau_c', alpha=1.0e-2, beta=1.0e-2)
    #Spatial clustering prior
        sigma = pm.Deterministic('sigma', np.sqrt(1/tau_c))
    #Co-efficents
        a0 = pm.Normal('a0', mu=0.0, tau=1.0e-3)
        a1 = pm.Normal('a1', mu=0.0, tau=1.0e-3)
        a2 = pm.Normal('a2', mu=0.0, tau=1.0e-3)
        a3 = pm.Normal('a3', mu=0.0, tau=1.0e-3)
    #Population Effect
        pop_eff= pm.Deterministic('pop_eff', tt.exp((-1*beta)/tt.exp(pop_den_all))) 
    #Regional Random Effects
        b_new = CAR('b_new', w=wmat, a=amat, tau=tau_c, shape=1425)
    #Mean/Expected Event Occurance Rate
        lamda = pm.Deterministic('lamda', tt.exp(a0 + a1*insta + a2*shear + a3*odd + b_new))
    #Actual Occurrence of Events
        Tlatent_new = pm.Poisson('Tlatent_new', mu=lamda, shape=1425)
    #Observed Event Counts
        Tobs_new = pm.Binomial('Tobs_new', n=Tlatent_new, p=pop_eff, shape=1425, observed=Tobs)
    #Initialization
        start1 = {'betaaux': 2., 'tau_c_log__': 2., 'a0': 1., 'a1': 1., 'a2': 1., 'a3': 1., 'Tlatent_new': Tlatent, 'b_new': b}
        step1 = pm.Metropolis([a0, a1, a2, a3, betaaux, beta, tau_c, b_new, Tlatent_new])

with model1:
     trace1 = merge_traces([pm.sample(draws=15000, step=[step1], tune=5000, start=start1, chain=i)
                        for i in range(1)])  

Моя модель работает, но вывод не выглядит правильным. В частности, "Tlatent_new" не обновляется с начального значения, которое я назначаю ему в "Tlatent". Если я не пытаюсь включить "pop_eff" в мою модель, то есть убрать строку "Tobs_new = pm.Binomial(" Tobs_new ", n=Tlatent_new, p=pop_eff, shape=1425, наблюдаемый =Tobs)" Я считаю, что " Tlatent_new'изменится от начальных значений, указанных в' Tlatent '. Я не могу понять, почему эта дополнительная строка помешает модели обновить 'Tlatent' или как я могу обойти это.

Любые предложения относительно того, что может быть проблемой, будут оценены.

Answer 1

Меня немного удивляет, что он застрял полностью, но Метрополис действительно не работает в более высоких измерениях. Я бы не ожидал, что это сработает хорошо.

Вы можете использовать NUTS, если найдете способ избавиться от дискретных переменных (наблюдаемые из них подходят). В этом примере, возможно, просто модель Tlatent как непрерывная переменная - Binomial работает с непрерывной n,

Пара других мелочей:

• Ранее параметры масштаба InverseGamma были очень распространены некоторое время назад, но, похоже, они могут показывать очень неудачное поведение. Если вы действительно хотите неинформативный предварительный, вы можете использовать Джеффри до использования log_sigma = pm.Flat(); sigma = tt.exp(log_sigma), Но в большинстве случаев гораздо лучше использовать HalfStudentT или HalfCauchy.
• Не нужно использовать tau, sd обычно приятнее читать.