##################################################################
# example of Monte Carlo "inverting" the cumulative distribution
#  -> Exercise: rewrite 'rfun()' making use of a binary search
#
# GdA 6/2/2015
##################################################################

# (unnormalized) function to sample
fun <- function(x) x^2

# some functions:  x^2
#                  (x-50)^2
#                  sqrt(x) 
#                  exp(-x/20)
#                  sin(x*(2*pi/100))^2

#-------------------------------------------------------------------
# different functions to invert the cumulative

# 1. using a very primitive method
rfunP <- function(x, Fx) {  
  ru <- runif(1)
  
  if( ru <= Fx[1] ) return(x[1])

  for (i in 2:length(x)) {
    if( (Fx[i-1] < ru) & (ru <= Fx[i]) ) return(x[i])
  }
}

# 2. using some peculiarities of R (faster than Method 1)
rfunR <- function(x, Fx) length(x) + 1 - length( Fx[runif(1) <= Fx] ) 

# 3. method based on binary search -> left as exercise
# ....
# ....
#-------------------------------------------------------------------

# select method to use
rfun <- rfunR
    
x  <- 1:100       # possible discrete values of the variable
fx <- fun(x)      # unnormalized f(x)
fx <- fx/sum(fx)  # normalization
Fx <- cumsum(fx)  # cumulative, F(x) 

print(round(fx, 4))   # print() is needed if you 'source' the file
print(round(Fx, 4))

# test on n events
n <- 100000
xr <- rep(0, n)
for (i in 1:n)  xr[i] <- rfun(x, Fx)

cat( sprintf("Mean: %.2f\n", mean(xr)) )
cat( sprintf("Std:  %.2f\n", sd(xr)) )

barplot(table(xr), col='cyan', main=body(fun))