# da 6 a 101 scatole ...
#
# G. D'Agostini, PED, nov. 2009
#---------------------------------------------------

n=101
prior.H       <- rep(1/n,n)        # prior
pW.dato.H     <- c(0:(n-1))/(n-1)  # likelihood W
pB.dato.H     <- 1 - pW.dato.H     # likelihood B

# Theorema di Bayes
bayes.nbox <- function(pObs.dato.H, p.H) {
  # pObs.dato.H e' il 'vettore' delle prob. condizionate
  # p.H         e' il vettore delle prior
  posteriors = pObs.dato.H * p.H / sum(pObs.dato.H * p.H)
  return(posteriors)
}

analizza.seq <- function(ris, p.H) {
  nr <- length(ris)
  if (nr == 1) ris <- as.vector(ris) 
  for (i in 1:nr) {
    if (ris[i] == 0 || ris[i] == 'B' || ris[i] == 'b') {
      p.H <- bayes.nbox(pB.dato.H, p.H)
    } else {
      p.H <- bayes.nbox(pW.dato.H, p.H)
    }
  }
  return(p.H)
}

ris <- c(1, 0, 0, 0, 1, 0, 0, 0, 0, 1)
p.H <- analizza.seq(ris, prior.H)
plot(p.H)

# ovviamente le probabilita' di H_i sono le stesse
# delle n proporzioni di palline bianche, ovvero di pW.dato.H
# che ne seguito possiamo chiamare semplicemente p:
p <- pW.dato.H 

# quindi
plot(p, p.H)

# memorizziamo questo risultato in
p.H.10 <- p.H
# e procediamo, ad es.
ris <- c(ris, c(0, 1, 0, 1, 0, 0, 0, 1, 0, 0) )  # aggiungiamo 10 risultati
p.H.20 <- analizza.seq(ris, prior.H)

ris <- c(ris, c(0, 0, 0, 0, 0, 0, 1, 0, 0, 0) )  # altri 10
p.H.30 <- analizza.seq(ris, prior.H)

ris <- c(ris, c(0, 0, 0, 1, 0, 1, 0, 0, 0, 0) )  # altri 10
p.H.40 <- analizza.seq(ris, prior.H)

ris <- c(ris, c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0) )  # altri 10
p.H.50 <- analizza.seq(ris, prior.H)

ris <- c(ris,  c(0, 0, 0, 0, 1, 0, 1, 0, 0, 0) )  # altri 10
p.H.60 <- analizza.seq(ris, prior.H)

ris <- c(ris, c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0) )  # altri 10
p.H.70 <- analizza.seq(ris, prior.H)

ris <- c(ris, c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0))  # altri 10
p.H.80 <- analizza.seq(ris, prior.H)

ris <- c(ris, c(0, 1, 0, 1, 0, 0, 0, 1, 0, 0))  # altri 10
p.H.90 <- analizza.seq(ris, prior.H)

ris <- c(ris, c(0, 0, 1, 0, 0, 0, 0, 0, 1, 1))  # altri 10
p.H.100 <- analizza.seq(ris, prior.H)

plot(p, p.H.100, col='blue', ylab='P(p)')
points(p, p.H.90, col='green')
points(p, p.H.80, col='orange')
points(p, p.H.70, col='red')
points(p, p.H.60, col='brown')
points(p, p.H.50, col='blue')
points(p, p.H.40, col='green')
points(p, p.H.30, col='orange')
points(p, p.H.20, col='red')
points(p, p.H.10, col='brown')
points(p, prior.H, col='black')

# chiaramente, 101 scatole danno l'idea di cosa
# succede per n->Inf
plot(p, p.H.100, col='blue', ty='l', ylab='f(p)')
points(p, p.H.90, col='green', ty='l', ylab='f(p)')
points(p, p.H.80, col='orange', ty='l', ylab='f(p)')
points(p, p.H.70, col='red', ty='l', ylab='f(p)')
points(p, p.H.60, col='brown', ty='l', ylab='f(p)')
points(p, p.H.50, col='blue', ty='l', ylab='f(p)')
points(p, p.H.40, col='green', ty='l')
points(p, p.H.30, col='orange', ty='l')
points(p, p.H.20, col='red', ty='l')
points(p, p.H.10, col='brown', ty='l')
points(p, prior.H, col='black', ty='l')

# essendo numeri oordinati, ha senso calcolare
# valore atteso e dev. standard di p.
# ad esempio dopo la 50-ma prova:
E.p <- sum( p * p.H.100 )  # e' esattamente la prob. di W in una prova futura
E.p

# confrontiamo questo valore con la frequenza relativa di 'successo'
f.W.100 <- length(ris[ris==1]) / length(ris)
f.W.100

# si confronti
E.p
# con
(length(ris[ris==1]) + 1) / (length(ris)+2)  # sorprendente?

#---------------------------------------------------
# a questo punto si puo' porre il nr di scatole a 1000 o 10000
# e osservare che il risultato e' molto stabile


# provare inoltre cosa succede se si osservano solo palline bianche
# o sole palline nere, ad es
ris = rep(0,10)
p.H <- analizza.seq( ris, prior.H)
plot(p, p.H, ylim=c(0,0.35),col='brown', ty='l', ylab='f(p)')
points(p, prior.H, ty='l')

p.H <- analizza.seq( c(ris,ris), prior.H)
points(p, p.H, ty='l',col='red')

p.H <- analizza.seq( c(ris,ris,ris), prior.H)
points(p, p.H, ty='l',col='green')

p.H <- analizza.seq( c(ris,ris,ris,ris), prior.H)
points(p, p.H, ty='l',col='blue')