#---------------------------------------------- # linear fit including inference on sigma and extrapolation # (same as done in OpenBUGS) # # GdA 27 marzo 2015 #---------------------------------------------- # function pausa pausa <- function() { cat ("\n >> Guarda il plot e dai enter per continuare\n") scan() } library(rjags) # list of 'data' data <- NULL m <- 2; c <- -1; sy <- 1.2 mux <- 3:20 muy <- m*mux + c data$n <- length(mux) data$x <- mux data$y <- muy + rnorm(data$n, 0, sy) data$x1 <- 30 plot(data$x, data$y, col='blue') fit <- lm(data$y ~ data$x) abline(fit, col='red') # list of inits inits <- list( m = 1 , c = 0 , tau = 0.25 , y1 = 20 ) # model model <- "fit_model.txt" jm <- jags.model(model, data, inits) update(jm, 100) # burn in to.sample <- c('m', 'c', 'sigma', 'muy1', 'y1') # variables of interest catena <- coda.samples(jm, to.sample, n.iter=20000) # sampling print(summary(catena)) pausa() # transform the chain in a dataframe catena.df <- as.data.frame( as.mcmc(catena) ) hist(catena.df$m, nc=100, col='blue') # etc. etc. pausa() plot(catena.df$m, ty='l', col='blue') pausa() plot(catena.df)