################################################################## # example of Monte Carlo reweighing uniformly extracted # possible values with the probability function # # GdA 9/2/2015 ################################################################## # (unnormalized) function to sample fun <- function(x) (x-50)^2 # some functions: x^2 # (x-50)^2 # sqrt(x) # exp(-x/20) # sin(x*(2*pi/100))^2 # select method to use rfunw <- function(x, fx, fxM) { accepted = FALSE while( !accepted ) { xr <- sample(x)[1] if ( runif(1)*fxM <= fx[x==xr] ) accepted = TRUE } return(xr) } x <- 1:100 # possible discrete values of the variable fx <- fun(x) # unnormalized f(x) fx <- fx/sum(fx) # normalization fxM <- max(fx) 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] <- rfunw(x, fx, fxM) E.X <- sum( x * fx ) sigma.X <- sqrt( sum( x^2 * fx ) - E.X^2 ) cat( sprintf("\nProbability distribution:\n") ) cat( sprintf(" Expected value: %.2f\n", E.X) ) cat( sprintf(" Standard deviation: %.2f\n\n", sigma.X) ) cat( sprintf("Simulated sample:\n") ) cat( sprintf(" Mean: %.2f\n", mean(xr)) ) cat( sprintf(" Std: %.2f\n", sd(xr)) ) barplot(table(xr), col='cyan', xlab='x', ylab='occurrences', main=body(fun))