0

我已经绘制了一个代码，其中说明了一些贝叶斯线性模型，我想弄清楚如何评估其后验分布下某个参数的任何概率，例如如何评估 P(-0.01<b_0<0.01|x)。我已经说明的模型如下（gama inversa 是 inverse gama）：

在此处输入图像描述

library(R2OpenBUGS)
library(invgamma)

mymodel = function(){
  for (i in 1:N) {
    Y[i] ~ dnorm(mu[i],tau)
    mu[i] <- b0+b1*x1[i]+b2*x2[i]+b3*x3[i]
  }

  b0 ~ dnorm(0,100)
  b1 ~ dnorm(0,100)
  b2 ~ dnorm(0,100)
  b3 ~ dnorm(0,100)
  tau ~ dnorm(0,100)

}

data=list(Y=c(250,        150,        128,        134,        110,        131,        98,        84,        147,        124,        128,        124,        147,        90,        96,        120,        102,        84,        86,        84,        134,        128,        102,        131,        84,        110,        72,        124,        132,        137,        110,        86,        81,        128,        124,        94,        74,        89), x1=c(        81.69,        103.84,        96.54,        95.15,        92.88,        99.13,        85.43,        90.49,        95.55,        83.39,        107.95,        92.41,        85.65,        87.89,        86.54,        85.22,        94.51,        80.80,        88.91,        90.59,        79.06,        95.50,        83.18,        93.55,        79.86,        106.25,        79.35,        86.67,        85.78,        94.96,        99.79,        88.00,        83.43,        94.81,        94.94,        89.40,        93.00,        93.59), x2=c(64.5,        73.3,        68.8,        65.0,        69.0,        64.5,        66.0,        66.3,        68.8,        64.5,        70.0,        69.0,        70.5,        66.0,        68.0,        68.5,        73.5,        66.3,        70.0,        76.5,        62.0,        68.0,        63.0,        72.0,        68.0,        77.0,        63.0,        66.5,        62.5,        67.0,        75.5,        69.0,        66.5,        66.5,        70.5,        64.5,        74.0,        75.5), x3=c(118,        143,        172,        147,        146,        138,        175,        134,        172,        118,        151,        155,        155,        146,        135,        127,        178,        136,        180,        186,        122,        132,        114,        171,        140,        187,        106,        159,        127,        191,        192,        181,        143,        153,        144,        139,        148,        179),N=38)
inits=function() {list(b0=1,b1=1,b2=1,b3=1,tau=1)}
exit=bugs(data=data,inits=inits, parameters.to.save = c("b0","b1","b2","b3"),
           model.file=mymodel,
           n.chains=1, n.iter=10000)
4

0 回答 0