Coverting OpenBUGS代碼PYMC3

我以前使用OpenBUGS/WinBUGS軟件工作執行我的貝葉斯統計，但已經決定搬過來在Python使用PYMC3包。所以我對pacakage相當陌生，仍然在學習如何充分利用它。我在將BUGS代碼轉換爲PYMC3時遇到一些困難。原來BUGS代碼如下：Coverting OpenBUGS代碼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，形狀= 1425，實測值= TOBS）' 我發現' Tlatent_new「將從'Ttent'中給出的初始值改變。我不明白爲什麼這個附加線路將防止模型更新「Tlatent」或如何，我可以解決這個問題。

任何建議，有什麼問題可能將不勝感激。

來源

2017-06-01 A.Cosmo

這是一個有點出人意料，我認爲它完全被卡住，但大都會真的不高維工作。我不希望這會在很好的工作。

你可以使用NUTS如果你找到一種方法來擺脫離散變量（觀測者的罰款）。在這個例子也許只是模擬Tlatent作爲一個連續變量 - Binomial連續n工作。

幾個其他的小東西：

的InverseGamma之前對尺度參數是很常見的，前一段時間，但似乎它可以顯示非常不幸的行爲。如果你真的想要一個無信息之前，你可以先使用log_sigma = pm.Flat(); sigma = tt.exp(log_sigma)的傑弗裏的使用。但在大多數情況下，使用HalfStudentT或HalfCauchy要好得多。
沒有必要使用tau，sd通常是更好的閱讀。

來源

2017-06-01 08:30:11 aseyboldt

Coverting OpenBUGS代碼PYMC3

回答

相關問題