Задайте условную линейную гауссову сеть с помощью rjags

Я изо всех сил пытаюсь определить условную линейную гауссовскую байесовскую сеть, используя rjags. (clg BN определяется непрерывным дочерним узлом (результатом), имеющим как непрерывную нормаль, так и дискретный родитель (предикторы))

Для приведенной ниже сети A дискретна, D и E непрерывны:

введите здесь описание изображения

Я полагаю, что для модели rjags я хочу, чтобы параметры узла 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"))

person user2957945    schedule 23.08.2017    source источник


Ответы (1)


Вам просто нужно использовать переменную 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. Таким образом, вы можете настроить значения и длины векторов (при условии совпадения длин catprobs, int, наклонов и знаков) без изменения кода JAGS. Например (используя runjags для удобства, хотя также возможно и с jags):

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

Мэтт

person Matt Denwood    schedule 23.08.2017
comment
большое спасибо, это именно то, на что я надеялся. - person user2957945; 24.08.2017