# inferenza parametrica di mu e sigma
# dato un campione di misure
#
# G. D'Agostini      PED, 26/11/09
#-------------------------------------------


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

prior.piatta.mu.sigma <- function(mu, sigma) {
  # mettiamo una prior piatta, con sigma >0
  if (sigma <= 0) return(0)
  else            return(1)
}

nn.pdf <- function(mu, sigma, x) {
  # posterior non normalizzata
  
  # check prior
  pr <- prior.piatta.mu.sigma(mu, sigma)
  if (pr == 0) return(0)
        
  n <- length(x)
  if (n < 1) return(0) 
  return( sigma^(-n) * exp( - sum( (mu - x)^2 ) / (2*sigma^2) ) * pr )
}

pdf.mu.sigma <- function(mu, sigma, x) {
  # pdf correttamente normalizzata
  nn.pdf(mu, sigma, x)/norma
}

# generiamo il campione
n          <- 10
mu.vero    <- 7
sigma.vera <- 2
x <- rnorm(n, mu.vero, sigma.vera)

# stimiamo mu e sigma da media e deviazione standard
mu.0       <- mean(x)
sigma.x.0  <- sd(x)*sqrt(n-1)/sqrt(n)
sigma.mu.0 <- sigma.x.0/sqrt(n)

#-------------------------------------------------------------------
# calcoliamo l'integrale di normalizzazione
# 1. in base ai valori che ci aspettiamo, definiamo una opportuna
# griglia 100x100
ng <- 100
sigma.max <- sigma.x.0 * (1 + 7./sqrt(2*n))         # (*)
sigma.min <- max(0, sigma.x.0 * (1 - 3./sqrt(2*n))) # (*)
dsigma <- (sigma.max - sigma.min)/ng
si <- seq(sigma.min + dsigma/2, sigma.max - dsigma/2, by=dsigma)
# (*) Nota: il range e' calcolato sapendo piu' o meno quello
#           che uno si aspetta; ma questo non e' un vero problema:
#           -> osservando il plot della pdf finale si controlla
#              se la griglia e' stata ben scelta; nel caso contrario
#              la si cambia e si fa rigirare lo script
#              (l'altra possibilitare e' di farne una piu' grande,
#               a scapito del tempo di calcolo)

mu.max    <- mu.0 + 5 * sigma.mu.0    
mu.min    <- mu.0 - 5 * sigma.mu.0    
dmu <- (mu.max-mu.min)/ng
mi  <- seq(mu.min + dmu/2, mu.max - dmu/2, by=dmu)

# 2. facciamo l'integrale sommando sui volumi dei parallelepipedini
area.cella <- dsigma*dmu
int <- 0
for (i in 1:ng) {
  for (j in 1:ng) {
    int <- int + nn.pdf(mi[i], si[j], x) * area.cella
  }
}
norma <- int 

# 3. a questo punto nn.pdf()/norma ci da' la pdf congiunta
# -> pdf.mu.sigma()
#------------------------------------------------------------

# marginalizziamo

int <- 0
pdf.mu    <- rep(0,ng)  # pdf
F.mu      <- rep(0,ng)  # cumulativa 
pdf.sigma <- rep(0,ng)
F.sigma   <- rep(0,ng)
for (i in 1:ng) {
  for (j in 1:ng) {
    int <- int + pdf.mu.sigma(mi[i], si[j], x) * area.cella   
    pdf.mu[i]    <- pdf.mu[i] + pdf.mu.sigma(mi[i], si[j], x)*dsigma
    pdf.sigma[i] <- pdf.sigma[i] + pdf.mu.sigma(mi[j], si[i], x)*dmu
  }
  if(i==1) {
    F.mu[i]    <- pdf.mu[i]*dmu
    F.sigma[i] <- pdf.sigma[i]*dsigma 
  } else {
    F.mu[i]    <- F.mu[i-1] + pdf.mu[i]*dmu
    F.sigma[i] <- F.sigma[i-1] + pdf.sigma[i]*dsigma
  }
}

#controlliamo che le pdf siano normalizzate:
sum( pdf.mu*dmu )
sum( pdf.sigma*dsigma )
# ovvero
F.mu[ng]
F.sigma[ng]

# plottiamo il risultato
plot(mi, pdf.mu, xlab='mu', ylab='f(mu)', ty='l')
pausa()
plot(mi, F.mu, xlab='mu', ylab='F(mu)', ty='l')
pausa()

plot(si, pdf.sigma,  xlab='sigma', ylab='f(sigma)', ty='l')
pausa()
plot(si, F.sigma,  xlab='sigma', ylab='F(sigma)', ty='l')

# calcoliamo valore atteso e incertezza standard delle due variabili:
E.mu    <- sum( mi * pdf.mu*dmu )
Var.mu  <- sum( mi^2 *pdf.mu*dmu ) - E.mu^2
std.mu  <- sqrt(Var.mu)
moda.mu <- mi[pdf.mu==max(pdf.mu)]

E.sigma    <- sum( si * pdf.sigma*dsigma )
Var.sigma  <- sum( si^2 * pdf.sigma*dsigma ) - E.sigma^2
std.sigma  <- sqrt(Var.sigma)
moda.sigma <- si[pdf.sigma==max(pdf.sigma)]

# si confrontino:
# valori veri
print("valori 'veri' di mu e sigma (vedi script)")
print(mu.vero)
print(sigma.vera)

# --> mu
print("risultati su mu")
print(E.mu)
print(mu.0)

print(std.mu)
print(sigma.mu.0)

print(moda.mu)

# --> sigma
print("risultati su sigma")
print(E.sigma)
print(sigma.x.0)

print(std.sigma)
print(sigma.x.0/sqrt(2*n))   #(*)

print(moda.sigma)
#----------------------------------------------------
# (*) per ora questo sigma.x.0/sqrt(2*n) e' misterioso,
#     ma vediamo che 'somiglia a std.sigma'

# -> provare cosa succede per n 'grande' (es. 100)
# -> provare a cambiare la prior

# -> eventualmente provare anche a cambiare la griglia   

#    ATTENZIONE: non esagerare con n grande, in quanto 
#    nn.pdf() puo' andare in underflow! (serve qualche trucco
#    numerico, ma la cosa non ci interessa troppo)
#    -> comunque, dopo un po' si capisce quando si possono usare
#       metodi approssimativi