#--------------------------------------------------------------
# inferenza di mu della normale con successiva previsione
#  - Variante in cui sigma.p (legata a osservazioni 'passate')
#    e sigma.f (legata alle osservazioni 'future')
#    possono differire  (ma entrambe le SIGMA sono ritenute NOTE,
#     in questo modello semplificato)
#  - Nel caso le osservazioni derivano da medie aritmetiche   
#    è sufficiente considerare la media come una misura equivalente, 
#    con sigma.p ridotta di sqrt(n.p), con 'n.p' il numero di valori
#     osservato nel passato.
#  - Si procede in modo analogo se si pensa a una media futura, 
#    tenendo conto che 'n.p' e 'n.f' possono differire
#
#  GdA, maggio 2025 
#---------------------------------------------------------

sigma.p = 1    # deviazione standard associata all'osservazione   
xp      = 5    # valore osservato

sigma.f = 2    # deviazione standard associata ali valori 
               # che si potranno osservare in un esperimento futuro  

E.mu  <- xp
sd.mu <-  sigma.p 
cat(sprintf("\n E(mu) = %0.2f,   std(mu) = %.2f\n", E.mu, sd.mu))

# limiti per gli istogrammi
mu.lim <- E.mu +  3*sqrt(sd.mu^2 + sigma.f^2)*c(-1,1)

N = 1000000
# valori di mu estratti secondo f(mu|E.mu,sd.mu) 
mur <- rnorm(N, E.mu, sd.mu)

# valori di xf che tengono conto di 'tutti i possibili 'mu'
xf    <- rnorm(N, mur, sigma.f)
E.xf  <- mean(xf)
sd.xf <- sd(xf)
cat(sprintf("\n E(x.f) = %0.2f,   std(x.f) = %.2f\n", E.xf, sd.xf))

par(mfrow= c(2,1) )
hist(mur, nc=100, col='cyan',  prob=TRUE,
     xlab='mu', ylab='f ( mu )', xlim=mu.lim)
m <- seq(E.mu-4*sd.xf, E.mu+4*sd.xf, len=101)
points(m, dnorm(m, E.mu, sd.mu), ty='l',  col='blue')

hist(xf, nc=100, col='cyan',  prob=TRUE,
     xlab='x.f', ylab='f ( x.f )',  xlim=mu.lim)
points(m, dnorm(m, E.mu, sqrt(sd.mu^2 + sigma.f^2)), ty='l',  col='blue')
     
par(mfrow= c(1,1) )