#----------------------------------------
# Biliardo di Bayes
#
# GdA PED, 2/12/2010
#-----------------------------------------

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

# plot biliardo
plot.biliardo <- function(x0, fmax, w, h) {
  dev.off()
  dev.new(width=w,height=h)
  plot(x0, xlim=c(0,1), ylim=c(-2.6,fmax), pch='.', xlab='p=x0/l', ylab='f(p)', main="Biliardo di Bayes")
  #ymin = -2.5
  #ymax = -0.5
  dy = ymax - ymin 
  lines(c(0,1), c(ymin, ymin),)
  lines(c(0,1), c(ymax, ymax))
  lines(c(0,0), c(ymin, ymax))
  lines(c(1,1), c(ymin, ymax))
  lines(c(x0,x0), c(ymin, ymax))
}

# plot palle
plot.palle <- function(x0, x, y) {
  x.l <- x[x<=x0]
  y.l <- y[x<=x0] * dy + ymin
  x.r <- x[x>x0]
  y.r <- y[x>x0]  * dy + ymin
  points(x.l, y.l, col='black', cex=2, pch='.')
  points(x.r, y.r, col='red', cex=2, pch='.')  
}
plot.beta <- function(xp, r,s, col) {
  points(xp, dbeta(xp, r, s), ty='l' , col=col)
}

ymin = -3.5
ymax = -0.5
dy = ymax - ymin

#n.evts <- 1000
#x0 <- runif(1)
#x <- runif(n.evts)
#y <- runif(n.evts)
#plot.biliardo(x0)
#plot.palle(x0, x, y)
#pausa()

n.evts <- 1000
x <- numeric()
y <- numeric()
x0 <- runif(1)
# np = as.integer( c(1,3,10,30,100,300,1000) )
np = seq(1,as.integer(sqrt(n.evts)),by=3)^2
np = c(np[np<n.evts], n.evts)
xp <- seq(0,1,by=0.0001)
r.stored <- numeric()
s.stored <- numeric()
plot(xp)
plot.biliardo(x0, 30, 16, 10)
for (i in 1:n.evts) {
  x[i] = runif(1)
  y[i] = runif(1)
  plot.palle(x0, x, y)
  r = length(x[x<=x0]) + 1
  s = length(x[x>x0]) + 1
  # cat(sprintf("(r,s) = (%d,%d)\n", r,s)) 
  plot.beta(xp, r, s, "blue")
  rnorm(5000) # come pausa
  if(i < n.evts) plot.beta(xp, r, s, "white")
  if(any(np == i) ){
    lr = length(r.stored) + 1
    r.stored[lr] = r
    s.stored[lr] = s
    cat(sprintf("(r,s) = (%d,%d)\n", r,s)) 
  }
  for (ir  in 1:length(r.stored)) {
    plot.beta(xp, r.stored[ir], s.stored[ir], "black")
  }
  #score = sprintf("(r,s) = (%d,%d)\n", r,s)
  #text(0.5, -3.6, score, cex=1.5)
  #text(0.5, -3.6, score, cex=1.5, col='white')
}