# inferenza parametrica di mu e sigma
# dato un campione di misure
#
# G. D'Agostini      PED, 26/11/09 (aggiornato 14/12/10)
#-------------------------------------------


# 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
}


# dati PED 9/12/2010
n          <- 10
mu.vero    <- 3
sigma.vera <- 1
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

# in alternativa si puo' generare un campione,
# modificando le righe seguenti
nuovo.campione = FALSE
if (nuovo.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 ng*ng
ng <- 200
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()
#------------------------------------------------------------

# punto di massima probabilita' nel piano (mu,sigma)
f.max = -1
p.max = c(NA, NA)
z <- matrix(rep(0,length(mi)*length(si)), length(mi),length(si) )
for (i in 1:length(mi)) {
  for (j in 1:length(si)) {
    z[i,j] <-  nn.pdf(mi[i], si[j], x) / norma
  }
}

for (i in 1:length(mi)) {
  for (j in 1:length(si)) {
    if (z[i,j] > f.max) {
      f.max = z[i,j]
      p.max = c(mi[i], si[j])
    }
  }
}
cat(sprintf("\nfmax=%f in (%f,%f)\n\n",f.max, p.max[1], p.max[2]))



# contour plot della pdf bidimensionale

filled.contour(mi, si, z,
               plot.title = title(main = "inferenza congiunta di mu e sigma",
               xlab = "mu", ylab = "sigma"),
               key.title = title(main="f(mu,sigma)"),)    

pausa()


par(bg="lightyellow", cex=1.2)
persp(mi,si,z,theta=30,phi=30,expand=0.5,
      col="lightblue",ltheta=90,shade=0.5,
      ticktype="detailed", xlab="mu", ylab="sigma",
      zlab="f(mu,sigma)", main="Inferenza congiunta")

pausa()

persp(mi,si,z,theta=-50,phi=25,expand=0.5,
      col="lightblue",ltheta=90,shade=0.5,
      ticktype="detailed", xlab="mu", ylab="sigma",
      zlab="f(mu,sigma)", main="Inferenza congiunta")

pausa()


cat(sprintf("\nI plot 3D che seguono (fatti con Open GL)\n"))
cat(sprintf("sono INTERATTIVI e vengono aperti su un'altra finestra (RGL):\n"))
cat(sprintf("\n  => ingrandire la finestra, quindi ruotare il plot con il mouse\n\n"))
library(rgl)
x3 <- y3 <- z3 <- rep(0, length(mi)*length(si))
k <- 0
for (i in 1:length(mi)){
  for (j in 1:length(si)){
    k <- k+1
    x3[k] <- mi[i]
    y3[k] <- si[j]
    z3[k] <- z[i,j]
  }
}
plot3d(x3,y3,z3, lwd=4,type="l",col=rainbow(1000),xlab='mu', ylab='sigma', zlab='f(mu,sigma)')

pausa()

decorate3d(x3,y3,z3, lwd=4,type="l",col=rainbow(1000),xlab='mu', ylab='sigma', zlab='f(mu,sigma)')

pausa()

# mettiamo i punti in ordine crescente di z3 in modo tale
# da avere bande di colore che cambiano con f(mu, sigma)
z3.ind <- sort.int(z3, index=TRUE)$ix
x3.s <- x3[z3.ind]
y3.s <- y3[z3.ind]
z3.s <- z3[z3.ind]
plot3d(x3.s, y3.s, z3.s,
       lwd=4,type="l",col=rainbow(length(x3)/2),
       xlab='mu', ylab='sigma', zlab='f(mu,sigma)')
# provare a giocare con i parametri!

pausa()

cat(sprintf("\nI rimanenti plot continuano sull'usuale finestra grafica di R\n\n"))

# 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)")
cat( sprintf("valori 'veri' (del generatore): mu = %f, sigma = %f\n\n",
     mu.vero,sigma.vera))  
#print(mu.vero)
#print(sigma.vera)

# --> mu
print("risultati su mu:")
cat( sprintf("  E(mu) = %f,  std(mu) = %f, moda(mu) = %f \n",
             E.mu, std.mu, moda.mu)) 
cat( sprintf(" [ media(x) = %f,  std(x)/sqrt(n) = %f]\n\n", mu.0, sigma.mu.0))

# --> sigma
print("risultati su sigma:")
cat( sprintf("  E(sigma) = %f,  std(sigma) = %f, moda(sigma) = %f \n",
             E.sigma, std.sigma, moda.sigma)) 
cat( sprintf(" [ std(x) = %f,  std(x)/sqrt(2*n) = %f]\n",
             sigma.x.0, sigma.x.0/sqrt(2*n)))


#----------------------------------------------------
# (*) 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

#-----------------------------------------------------------------
# inoltre: provare a fare l'inferenza su mu assumendo
# un certo valore valore di sigma (ad es. sigma=1 usato per
# generare il campione)
# -> "inferenza condizionata"
#    -> concettualmente facile:
#       -> fissare un certo valore di sigma, sia sigma0,
#          nella f(mu,sigma|dati),
#          che a questo punto diventa nn.f(mu|sigma0,dati):
#          -> quindi rimane soltanto da rinormalizzare nn.f()
#              per ottenere f(mu|sigma0,dati)
#-----------------------------------------------------------------