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