#----------------------------------------------------------
# Semplice Gibbs sampler per inferire mu e sigma
# da un campione normale di dati, con prior uniformi
# (con distribuzione predittiva di ulteriore osservazione)
#
# Caso semplice, con prior piatte
#
# GdA 24/04/09
#-----------------------------------------------------------
 
# teoria:
# f(mu,sigma|{x}) \propto f({x}|mu,sigma)
#                 \propto sigma^(-n)*exp[-n/(2*sigma^2) * ((mu-x.m)^2+s^2)]
#                 \propto tau^(n/2)*exp[-n*tau/2 * ((mu-x.m)^2+s^2)]     
#
#  con: tau = 1/sigma^2
#       x.m: media aritmetica delle osservazioni
#       s^2: media dei quadrati degli scarti
#  [si ricorda la relazione: \Sum_i (x_i-mu)^2 = n*((mu-x.m)^2 + s^2) ]
#
# ne segue che:
#
# f(sigma|{x},mu)  \propto  tau^(n/2) * exp[-tau * (n/2)*((mu-x.m)^2+s^2) ] 
#
#                  -> gamma con parametri (n/2)+1 e (n/2)*((mu-x.m)^2+s^2)
#
# f(mu|{x},sigma) \propto  exp[-tau* n/2 * ((mu-x.m)^2 ]
#                  -> normale con media x.m e sigma=1/sqrt(tau*n)
#
# inoltre, per quanto riguarda la distribuzione predittiva:
#
# f(y|mu,sigma) -> normale com mu e sigma
 
# genera dati
 
mu     = 10
sigma  = 2
n.data = 1000    # interessante provare con n.data molto piccolo
 
x <- rnorm(n.data, mu, sigma)
 
# calcola summaries:
x.m = mean(x)
s  = sqrt(mean(x^2) - x.m^2)
n  = length(x) 
 
# parametri mcmc
n.sample = 10000
mu0    = 0
sigma0 = 1
 
 
mu    <- numeric()
sigma <- numeric()
tau   <- numeric()
y     <- numeric()
 
# Gibbs sampler
mu[1]    <- mu0
sigma[1] <- sigma0
tau[1]   <- 1/sigma0^2 
y[1]     <- rnorm(1, mu[1], sigma[1])
 
for (i in 2:n.sample) {
  mu[i]  <- rnorm(1, x.m, sigma[i-1]/sqrt(n))
  tau[i] <- rgamma(1, n/2 + 1,
                      n/2*((mu[i]-x.m)^2 + s^2) )
  sigma[i] <- 1/sqrt(tau[i])
  y[i]   <- rnorm(1, mu[i], sigma[i])
}
 
#---------------------------------------------------------
# Comandi per analizzare la catena
# 
# (attenzione: catena risente del fatto che il punto
# iniziale e' 'scelto male': -> ma converge
# dopo un paio di campionamenti.
plot(mu[1:20])
# Per evitare problemi
# facciamo i summary eliminando i primissimi punti)
 
mean(mu[10:n.sample])
sd(mu[10:n.sample])
# confrontiamo con
mean(x)
sd(x)/sqrt(n)
 
mean(sigma[10:n.sample])
sd(sigma[10:n.sample])
 
mean(y[10:n.sample])
sd(y[10:n.sample])
 
hist(mu[10:n.sample],nc=50)
hist(y[10:n.sample],nc=50)
hist(sigma[10:n.sample],nc=50)