#####################################################################
## Rejection sampling: Example 
## Goal Y~Beta(6,3)
## Proposal X ~ U(0,1)
## c = 3

## Author: Andrea Riebler <andrea.riebler@math.ntnu.no> (16.01.2015)
#####################################################################


## plot the target distribution together with the proposal
## distribution multiplied by c
c = 5
par(mfrow=c(1,2), lwd=2)
curve(dbeta(x,6,3), from=0, to=1, ylim=c(0, c+1), ylab="", col="blue")
abline(h=c, col="red")
legend("top", c("f(x) = Beta(6, 3)", "c*g(x)=c*U(0,1)"), 
  col=c("blue", "red"), lty=1, lwd=2)

finished = 0
while(finished == 0){  
  ## step-by-step illustration  
  ## 1.) sample from the proposal distribution
  sample.x = runif(1, 0, 1)

  ## 2.) Compute the acceptance probability alpha
  alpha = dbeta(sample.x, 6, 3)/(dunif(sample.x, 0, 1) * c)

  ## 3.) Generate a helping variable U from a uniform(0,1)
  u = runif(1, 0, 1)
  points(sample.x, dunif(sample.x, 0, 1) * c * u)
  Sys.sleep(1)
  ## 4.) Check if we accept it
  if(u <= alpha){
    y = sample.x
    col = "blue"
    points(sample.x, dunif(sample.x, 0, 1) * c * u, col=col)
    finished = 1
  } else {
    col = "red"
    points(sample.x, dunif(sample.x, 0, 1) * c * u, col=col)
  }
  Sys.sleep(2)
}    
  
  
## let us generate n realisations from the proposal distribution
n=10000
sample.x = runif(n, 0, 1)

## generate n realisations from a uniform to decide which of
## n proposed samples we keep, i.e. regard as realisations from
## Beta(6,3)
U = runif(n, 0, 1)

## compute the acceptance criterion
accept = (dunif(sample.x, 0, 1)*c*U) <= (dbeta(sample.x, 6, 3))

## check how many do we accept
table(accept)

## mark the points (sample.x, U*c) in the plot with 
## colors according whether they are accepted or not
points(sample.x[accept], U[accept]*c, col="blue", pch=20)
points(sample.x[!accept], U[!accept]*c, col="red", pch=20)

## check the distribution of all accepted points by making
## a histogram and compare to the target distribution.
library(MASS)
truehist(sample.x[accept], ylim=c(0,c+1), col="white")
curve(dbeta(x,6,3), from=0, to=1, add=TRUE)

## illustrate the distribution of all samples
par(mfrow=c(1,3), lwd=2)
truehist(sample.x, ylim=c(0,c+1), col="white", main="All generated samples")
truehist(sample.x[accept], ylim=c(0,c+1), col="white", main="Accepted samples")
truehist(sample.x[!accept], ylim=c(0,c+1), col="white", main="Rejected samples")