每个后验密度用计算机画出图形,使用R2WinBUGS拟合后怎么绘制后验概率密度

library(R2WinBUGS)

library(lattice)

library(coda)

dat

attach(dat)

setwd("c:/")

sink("lme.model2.txt")

cat("

model {

# Priors

for (i in 1:6){

alpha[i] ~ dnorm(mu.int, tau.int)        # Random intercepts

beta[i] ~ dnorm(mu.slope, tau.slope)# Random slopes

}

mu.int ~ dnorm(0, 0.001)                # Mean hyperparameter for random intercepts

tau.int

sigma.int ~ dunif(0, 100)                # SD hyperparameter for random intercepts

mu.slope ~ dnorm(0, 0.001)                # Mean hyperparameter for random slopes

tau.slope

sigma.slope ~ dunif(0, 100)                # SD hyperparameter for slopes

tau

sigma ~ dunif(0, 100)                        # Residual standard deviation

# Likelihood

for (i in 1:n) {

mass[i] ~ dnorm(mu[i], tau)

mu[i]

}

}

",fill=TRUE)

sink()

# Bundle data

win.data

# Inits function

inits

mu.int = rnorm(1, 0, 1), sigma.int = rlnorm(1), mu.slope = rnorm(1, 0, 1),

sigma.slope = rlnorm(1), sigma = rlnorm(1))}

# Parameters to estimate

parameters

ni

nb

nt

nc

# Start Gibbs sampling

out1

n.chains=nc, n.burnin=nb, n.iter=ni, debug = TRUE)