1

我以前一直在使用 OpenBUGS/WinBUGS 来执行我的贝叶斯统计，但我决定转而使用 Python 中的 PYMC3 包。所以我对 pacakage 还很陌生，还在学习如何充分使用它。我在将 BUGS 代码转换为 PYMC3 时遇到了一些困难。原BUGS代码如下：

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,observed=Tobs)' 我发现' Tlatent_new' 将从'Tlatent' 中给出的初始值改变。我不明白为什么这条附加线会阻止模型更新“潜在”或如何解决这个问题。

任何关于问题可能是什么的建议将不胜感激。

4

1 回答 1

1

我有点惊讶它完全卡住了，但 Metropolis 在更高维度上确实不起作用。我不希望这能很好地工作。

如果您找到摆脱离散变量的方法（观察到的变量很好），您可以使用 NUTS。在此示例中，可能只是将模型建模Tlatent为连续变量——Binomial与 Continuous 一起使用n。

其他几件小事：

  • 尺度参数的 InverseGamma 先验在前段时间很常见，但它似乎可以表现出非常不幸的行为。如果您真的想要一个无信息的先验，您可以使用 Jeffrey 的先验使用log_sigma = pm.Flat(); sigma = tt.exp(log_sigma). 但在大多数情况下，最好使用 HalfStudentT 或 HalfCauchy。
  • 无需使用tau，sd通常更易于阅读。
于 2017-06-01T08:30:11.567 回答