##---------------------------------------------------- # Gibbs sampling per inferire mu e sigma e predire # l'eventuale risultato successivo # # GdA 15 marzo 2012 ##---------------------------------------------------- source("pausa.R") ricampiona = TRUE verb=!TRUE isto=TRUE n.gibbs <- 200 n.sample = 10 # true parameters mu0 <- 1 sigma0 <- 1 # simulated data if( ricampiona ) { x <- rnorm(n.sample, mu0, sigma0) } ## check: cat( sprintf("%d valori: media e sigma: %f, %f\n", n.sample, mean(x), sd(x) ) ) # hist (x, nc=20) #----------------------------------------------------- # inference with Gibbs sampler # # 1) PRIORS ----------------------------------------------------- # very broad normal prior for mu is assumed mu.p <- 0 sigma.p <- 1000 # very broad Gamma prior for inverse of variance ('tau' in BUGS), # (remember: E[X]=c/r, var[X]=c/r^2 and that, if c=1, # the exponential is recovered) c.p <- 1 # -> exponential r.p <- 0.0001 # -> very flat! # 2) inits # initial value of mu and sigma mu.init <- 2 # 10 sigma.init <- 2 # 10 # 3) campionamento # vettori per la storia, inizializzatia zero h.mu <- h.sigma <- h.y <- rep(0, n.gibbs+1) mu <- mu.init sigma <- sigma.init h.mu[1] <- mu h.sigma[1] <- sigma h.y[1] <- h.mu[1] nx <- length(x) sx <- sum(x) xm <- mean(x) s2 <- var(x)*(nx-1)/nx # att! -> "(n-1)" non c'entra niente! for (i in 2:(n.gibbs+1)) { # sample mu from a gaussian (given the old sigma) E.mu <- ( sigma.p^2 * sx + sigma^2 * mu.p ) / (nx*sigma.p^2 + sigma^2) S.mu <- sqrt ( sigma^2 * sigma.p^2 / ( nx*sigma.p^2 + sigma^2 ) ) mu <- rnorm(1, E.mu, S.mu) # sample tau from a Gamma (given the old mu) c <- c.p + nx/2 # r <- r.p + sum((x-mu)^2)/2 r <- r.p + nx/2*(s2 + (xm-mu)^2) tau <- rgamma(1, c, r) sigma <- 1/sqrt(tau) h.mu[i] <- mu h.sigma[i] <- sigma h.y[i] <- rnorm(1, h.mu[i], h.sigma[i]) if(verb) { cat (sprintf ("E.mu %f, S.mu %f, mu %f; c %f, r %f, tau %f, sig %f, y %f\n", E.mu, S.mu, mu, c, r, tau, sigma, h.y[i] )) } } if(isto) { # plot(h.mu, h.sigma, cex=0.2) # pausa() plot(NULL, xlab='mu', ylab='sigma', xlim=c(min(h.mu),max(h.mu)), ylim=c(min(h.sigma),max(h.sigma))) n.p <- n.gibbs points(h.mu[1], h.sigma[1],cex=0.6, col='blue',pch=19) for(i in 1:(n.p-1)) { lines(h.mu[i:(i+1)], h.sigma[i:(i+1)], col='red') points(h.mu[i], h.sigma[i],cex=0.2) } points(h.mu[n.p], h.sigma[n.p],cex=0.6, col='green',pch=19) }