1

以前、OpenBUGS/WinBUGS を使用してベイジアン統計を実行していましたが、Python で PYMC3 パッケージを使用することにしました。そのため、私はパッケージにかなり慣れておらず、まだ完全に使用する方法を学んでいます. 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' で指定された初期値から変更されます。この追加の行によってモデルが「Tlatent」を更新できない理由、またはこれを回避する方法がわかりません。

問題が何であるかについての提案をいただければ幸いです。

4

1 に答える 1

1

完全に行き詰まってしまうのはちょっと驚きですが、Metropolis は実際には高次元では機能しません。これがうまく機能するとは思わないでしょう。

離散変数を取り除く方法が見つかった場合は、NUTS を使用できます (観測された変数は問題ありません)。この例Tlatentでは、連続変数としてモデル化するだけかもしれません –Binomial連続で動作しnます。

その他のいくつかの小さなこと：

  • スケール パラメーターの InverseGamma 事前設定は、以前は非常に一般的でしたが、非常に残念な動作を示す可能性があるようです。本当に有益でない事前確率が必要な場合は、ジェフリーの以前の使用を使用できますlog_sigma = pm.Flat(); sigma = tt.exp(log_sigma)。しかし、ほとんどの場合、HalfStudentT または HalfCauchy を使用する方がはるかに優れています。
  • tauを使用する必要はありませんsd。通常は読みやすいです。
于 2017-06-01T08:30:11.567 に答える