############################################################
## Illustration Monte-Carlo integration
############################################################

## Binomial distribution
###########################

n = 10000

x = rbinom(n, size=20, p=0.3)

## get the mean
(my.mean = mean(x))

## the mean is given by size*p
(20*0.3)

## get the variance
(my.var = mean(x^2) - mean(x)^2)

## compare with theoretical result
(20*0.3*0.7)

## beta distribution
################################

n = 1000

alpha=2
beta=4

x = rbeta(n, shape1=alpha, shape2=beta)
library(MASS)
truehist(x)
curve(dbeta(x, shape1=alpha, shape2=beta), from=0, to=1, add=T)

## get the mean
(my.mean = mean(x))

## the mean is given by alpha/(alpha+beta)
(alpha/(alpha+beta))

## get the variance alpha*beta/((alpha+beta)^2*(alpha+beta+1))
(my.var = mean(x^2) - mean(x)^2)

## compare with theoretical result
(alpha*beta/((alpha+beta)^2*(alpha+beta+1)))


############################################################
## use importance sampling and sample from a uniform
############################################################

x.help = runif(n)

## compute the importance weights
my.weight = dbeta(x.help, shape1=alpha, shape2=beta)/dunif(x.help)

## estimate the mean of the Beta(alpha, beta) distribution
(sum(my.weight*x.help)/sum(my.weight))

## estimate the variance of the Beta(alpha, beta) distribution
(sum(my.weight*x.help^2)/sum(my.weight) - (sum(my.weight*x.help)/sum(my.weight))^2)