#---------------------------------------
# primo esempio di inferenza gaussiana,
# assumendo sigma nota
#
# GdA PED, 13/12/2010
#---------------------------------------

# function pausa
pausa <- function() {
  cat ("\n >> Guarda il plot e dai enter per continuare\n")
  scan()
}

#---------------------------------------------------

prior.piatta.mu <- function(mu) {
  # mettiamo una prior piatta ('impropria')
  return(1)
}

lik.norm.sigma.nota <- function(mu, x, sigma) {
  lik <- 1
  for (i in 1:length(x)) {
    lik <- lik * dnorm(x[i], mu, sigma)
  }
  return(lik)
}

nn.pdf <- function(mu, x, sigma) {
  # posterior non normalizzata
  return( lik.norm.sigma.nota(mu, x, sigma) * prior.piatta.mu(mu) ) 
}

#----------------------------------------------
# dati PED 9/12/2010
x.ped = c(2.5, 1.7, 1.9, 2.8, 2.1, 2.5, 1.1, 2.8, 2.5, 2.1)
x = x.ped

mu.pl <- seq(mean(x)-3*sd(x), mean(x)+3*sd(x), len=1000)
plot(mu.pl, nn.pdf(mu.pl, x, 1),ty='l')

for (i in 1:length(x)){
  pausa()
  titolo = sprintf("inferenza con %d valore(i)",i)
  plot(mu.pl, nn.pdf(mu.pl, x[1:i], 1),ty='l', main=titolo, xlab='mu', ylab='~ f(mu)')
}

# rifacciamolo, simulando al volo n dati,
x <- rnorm(101, 2, 1)
mu.pl <- seq(mean(x)-3*sd(x), mean(x)+3*sd(x), len=3000)
for (i in seq(1,length(x),by=5) ){
  pausa()
  titolo = sprintf("NUOVO CAMPIONE: inferenza con %d valore(i)",i)
  plot(mu.pl, nn.pdf(mu.pl, x[1:i], 1),ty='l', main=titolo, xlab='mu', ylab='~ f(mu)')
}


# riprendiamo i dati originali e facciamo anche la normalizzazione
x <- x.ped

# integrale fatto numericamente (metodo "rozzo ma efficace")
mu.i = seq(mean(x)-3*sd(x), mean(x)+3*sd(x), len=10000)
norma = 0
dmu = mu.i[2]-mu.i[1]
dmu.h = dmu/2
for (i in 1:(length(mu.i)-1)) {
  norma <- norma +  nn.pdf(mu.i[i]+dmu.h, x, 1) * dmu
}

E.mu = 0
for (i in 1:(length(mu.i)-1)) {
  E.mu <- E.mu +  (mu.i[i]+dmu.h) * nn.pdf(mu.i[i]+dmu.h, x, 1)/norma * dmu
}

E.mu2 = 0
for (i in 1:(length(mu.i)-1)) {
  E.mu2 <- E.mu2 +  (mu.i[i]+dmu.h)^2 * nn.pdf(mu.i[i]+dmu.h, x, 1)/norma * dmu
}

sigma.mu = sqrt(E.mu2-E.mu^2)