2017-08-23 69 views
2

我正在努力使用rjags來定義條件線性高斯貝葉斯網絡。 使用rjags定義條件線性高斯網絡

對於下面的淨(A CLG BN是通過具有兩個連續的正常和離散父(預測)一個連續的子節點(結果)所定義的),A是離散的,d和E連續:

enter image description here

對於rjags模型,我supose我要的是節點E的參數進行數值節點A上定義需要:僞代碼

model { 
    A ~ dcat(c(0.0948, 0.9052)) 
    D ~ dnorm(11.87054, 1/1.503111^2) 

    if A==a then E ~ dnorm(6.558366 + 1.180965*D, 1/2.960002^2) 
    if A==b then E ~ dnorm(3.370021 + 1.532289*D, 1/6.554402^2) 
} 

我可以通過使用下面的代碼得到一些工作,但它會很快與更多的預測器和分類級別混淆。

library(rjags) 

model <- textConnection("model { 
    A ~ dcat(c(0.0948, 0.9052)) 
    D ~ dnorm(11.87054, 1/1.503111^2) 

    int = 6.558366 - (A==2)*(6.558366 - 3.370021) 
    slope = 1.180965 - (A==2)*(1.180965 - 1.532289) 
    sig = 2.960002 - (A==2)*(2.960002 - 6.554402) 

    E ~ dnorm(int + slope*D, 1/sig^2) 
}") 

jg <- jags.model(model, n.adapt = 1000 

我的問題:我該如何定義簡潔這種模式嗎?

數據從

來到
library(bnlearn) 
net = model2network("[A][D][E|A:D]") 
ft = bn.fit(net, clgaussian.test[c("A", "D", "E")]) 

coef(ft) 
structure(list(A = structure(c(0.0948, 0.9052), class = "table", .Dim = 2L, .Dimnames = list(
    c("a", "b"))), D = structure(11.8705363469396, .Names = "(Intercept)"), 
    E = structure(c(6.55836552742708, 1.18096500477159, 3.37002124328838, 
    1.53228891423418), .Dim = c(2L, 2L), .Dimnames = list(c("(Intercept)", 
    "D"), c("0", "1")))), .Names = c("A", "D", "E")) 

sigma(ft) 
structure(list(A = NA, D = 1.50311121682603, E = structure(c(2.96000206596326, 
6.55440224877698), .Names = c("0", "1"))), .Names = c("A", "D", 
"E")) 

回答

2

你只需要使用您的變量A作爲索引參數:

library('rjags') 

model <- " 
model { 
    A ~ dcat(c(0.0948, 0.9052)) 
    D ~ dnorm(11.87054, 1/1.503111^2) 

    ints <- c(6.558366, 3.370021) 
    int <- ints[A] 
    slopes <- c(1.180965, 1.532289) 
    slope <- slopes[A] 
    sigs <- c(2.960002, 6.554402) 
    sig <- sigs[A] 

    E ~ dnorm(int + slope*D, 1/sig^2) 
} 
" 

jg <- jags.model(textConnection(model), n.adapt = 1000) 

順便說一句,因爲你有很多在固定數量的模型,在R中定義它們然後將它們作爲數據傳遞給JAGS可能更有意義。這樣,您可以在不修改JAGS代碼的情況下調整矢量的值和長度(只要catprobs,ints,slope和sigs匹配)。例如(使用runjags爲方便起見,雖然也可以用尖齒):

library("runjags") 

model <- " 
model { 
    A ~ dcat(catprobs) 
    D ~ dnorm(Dmu, Dprec) 

    int <- ints[A] 
    slope <- slopes[A] 
    sig <- sigs[A] 

    E ~ dnorm(int + slope*D, 1/sig^2) 

    #data# catprobs, Dmu, Dprec, ints, slopes, sigs 
    #monitor# A, D, E 
} 
" 

catprobs <- c(0.0948, 0.9052) 
Dmu <- 11.87054 
Dprec <- 1/1.503111^2 
ints <- c(6.558366, 3.370021) 
slopes <- c(1.180965, 1.532289) 
sigs <- c(2.960002, 6.554402) 

results <- run.jags(model) 
results 

馬特

+0

非常感謝你,那正是我所期待的。 – user2957945