OpenBUGSコードをPYMC3にカバーする

Question

以前はOpenBUGS/WinBUGSを使用してベイジアン統計を行ってきましたが、PythonでPYMC3パッケージを使用することにしました。だから私はかなり新しい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、observed = 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は通常、読みやすいです。