# Funzione di probabilita'  
f.theta.prior  <- function(x) {
   y <- numeric()
   for (i in seq(along.with=x)) {
     if (x[i] < 1. || x[i] > 10.)  y[i] = 0
     else  y[i] <- 1./(log(10)*x[i])
   }
   return(y)
}
 
# guardiamo forma della pdf
curve(f.theta.prior(x),0,11)
# controlliamo la normalizzazione:
integrate( f.theta.prior, 0, 11)
 
# Inversione della cumulativa 
r.theta.prior  <- function(n) {
   l10 <- log(10)
   y <- numeric()
   x <- runif(n)
   for (i in seq(1:n)) {
     y[i] <- exp(l10*x[i])
   }
   return(y)
}
 
# hit/miss
r.theta.prior.hm  <- function(n) {
   fmax <- f.theta.prior(1)
   i <- 0
   r <- numeric()
   while (i < n) {
     x <- runif(1)*9 + 1
     y <- runif(1) * fmax
     if ( y <= f.theta.prior(x) ) {
        i    <- i + 1
        r[i] <- x 
     } 
   }
   return(r)
}
 
# usiamo primo metodo
theta.r <- r.theta.prior(10000)
hist(theta.r, nclass=100, prob=TRUE)
curve(f.theta.prior(x), 0, 11, n=200, add=TRUE)
 
# secondo metodo: OK, solo piu' lento
theta.r <- r.theta.prior.hm(10000)
hist(theta.r, nclass=100, prob=TRUE)
curve(f.theta.prior(x), 0, 11, n=200, add=TRUE)
 
#-------------------------------------------------------------
# Approfittiamo per fare un primo integrale con
# metodi di MC
#
# Integrale di  f.theta.prior fra 2 e 5
 
# 1) usando integrate()
th.min <- 2
th.max <- 5
integrate(f.theta.prior, th.min, th.max)   # -> 0.39794
 
# 2) integrazione con metodo di MC
x  <- runif(10000, th.min, th.max)
fx <- f.theta.prior(x)
sum(fx) * (th.max-th.min)/length(x)        # -> 0.3965 in un run
 
# Il primo metodo e' molto piu' accurato del secondo,
# ma il secondo puo' essere comodo in problemi multidimensionali