#------------------------------------------------------
#  Example of simulation based on hit/miss method
#
#  GdA 13/2/2015
#-------------------------------------------------

# function pausa  (useful to use this file as a script "sourcing" it)
pausa <- function() {
  cat ("\n >> Guarda il plot e dai enter per continuare\n")
  scan()
}

# unnormalized function to sample from
fun <- function(x) (sin(x)*x)^2    
xmin <- 0
xmax <- pi

# plot fun
curve(fun, xmin, xmax)
pausa()

n=100000

xu <- xmin + runif(n)*(xmax-xmin)   # uniform in the range

fxM <- max( fun(xu) )                # find maximum

yu <- runif(n)*fxM*1.1              # 1.1 is an (exagerated!) security factor

cond <- yu <= fun(xu)               # condition "hit"

plot(xu, yu, cex=0.1)               # plot all points

points(xu[cond], yu[cond], cex=0.1, col='red')  # only points satisfying 'cond'
# uncomment to also show fun()  
# curve(fun, xmin, xmax, add=TRUE, col='blue', lwd=2)

pausa()

xr <- xu[cond]                      # sample of interest
hist(xr, nc=200, col='red')

# Byproduct:

accept <- length(xr)/n
cat( sprintf("\n\n accepted points %.4f \n", accept) ) 

A <- accept*(xmax-xmin)*fxM*1.1    # Att! do not forget 1.1 !!

cat( sprintf(" Area below fun(x):  %.3f \n", A) )

A.int <- integrate(fun, xmin, xmax)
cat( sprintf(" integrate(fun, xmin, xmax) : %.6f \n\n", A.int[[1]]) )

exact.int <- pi*(2*pi^2-3)/12
cat( sprintf("[exact integral = %.6f] \n\n", exact.int) )

# to understand the object "A.int" (from console)
str(A.int)