## Estimate P(X > 1) where X follows N(0,1)

 1- pnorm(1)

 ## Using Monte Carlo Method
 
 ## P(X > 5) = integral I(x> 5) f(x) dx where f is the density of N(0,1)
 
 ## Draw Samples from N(0,1) and then take the average of the indicators, which 
 ## is the proportion of data larger than 5. 
 
 n=10000;
 x <- rnorm(n);
 sum(x > 5)/n
 
 ## Importance sampling method
 ## First draw samples from g, call them x_1, ..., x_n
 ## Then approximate the integral by average of f(x_i) h(x_i)/ g(x_i)
 ## Choose g to be normal with mean 5
 
 n=10000;

 
 ## f = N(0,1), g= N(5,1), h = I(x > 5)
 
 x <- rnorm(n,5,1);
 
 y <- (dnorm(x)/dnorm(x,5))*I(x>5)
 
 sum(y)/n
 
 ## f = N(0,1), g= N(6,1), h = I(x > 5)
 
 x <- rnorm(n,6,1);
 
 y <- (dnorm(x)/dnorm(x,6))*I(x>5)
 
 sum(y)/n
 
 ## f = N(0,1), g= N(5,4), h = I(x > 5)
 
 x <- rnorm(n,5,sd=2);
 
 y <- (dnorm(x)/dnorm(x,5,sd=2))*I(x>5)
 
 sum(y)/n
 
 ### REjection sampling
 # density function for f(x)
 densFun <- function(x) {return(sin(pi*x)^2)}
 x <-seq(0, 1, 10^-2)
 
 ## Plot densFun and show the uniform density proposal as envelope
 
 plot(x, densFun(x), 
      xlab = "x", ylab = "pdf(x)",
      type = "l", lty = 1, lwd = 2, col = "blue",
      main = "Densities comparison")
 lines(x, dunif(x, 0, 1), lwd = 2, col="red")
 legend('topright',
        c(expression(sin(pi*x)^2), "uniform(0,1)"), 
        lty = 1, lwd = 2, col = c("blue", "red"))
 
 ## Plot densFun and show the Beta density as envelope
 plot(x, densFun(x), 
      xlab = "x", ylab = "pdf(x)",
      type = "l", lty = 1, lwd = 2, col = "blue",
      main = "Densities comparison")
 lines(x, (1/1.5)*dbeta(x, 2, 2), lwd = 2, col="red")
 legend('topright',
        c(expression(sin(pi*x)^2), "uniform(0,1)"), 
        lty = 1, lwd = 2, col = c("blue", "red"))
 
 
 #c=1, q(x) = 1$
 
 numSim=10^4
 samples = NULL
 for (i in 1:numSim) {
   # get a uniform proposal
   proposal <- runif(1) 
   # calculate the ratio
   densRat <- densFun(proposal)/dunif(proposal)  
   #accept the sample with p=densRat
   if ( runif(1) < densRat ){ 
     #fill our vector with accepted samples
     samples <- c(samples, proposal) 
   }
 }
 
 #c= 1/1.5; q(x) = Beta (2,2)
 c=1/1.5;
 numSim=10^4
 samples = NULL
 for (i in 1:numSim) {
   # get a uniform proposal
   proposal <- rbeta(1,2,2) 
   # calculate the ratio
   densRat <- densFun(proposal)/(c*dbeta(proposal,2,2))  
   #accept the sample with p=densRat
   if ( runif(1) < densRat ){ 
     #fill our vector with accepted samples
     samples <- c(samples, proposal) 
   }
 }
 
 hist(samples, freq=FALSE) #construct density hist
 print(paste("Acceptance Ratio:", length(samples)/numSim))
 
 ### Write the rejection sampler as a function of the target and envelope
 
 sim_fun <- function(f, envelope = "unif", par1 = 0, 
                     par2 = 1, n = 10^2, plot = TRUE){
   
   r_envelope <- match.fun(paste0("r", envelope))
   d_envelope <- match.fun(paste0("d", envelope))
   proposal <- r_envelope(n, par1, par2)
   density_ratio <- f(proposal) / d_envelope(proposal, par1, par2)
   samples <- proposal[runif(n) < density_ratio]
   acceptance_ratio <- length(samples) / n
   if (plot) {
     hist(samples, probability = TRUE, 
          main = paste0("Histogram of ", n, " samples from ", 
                        envelope, "(", par1, ",", par2, ").\n 
                       Acceptance ratio: ", 
                        round(acceptance_ratio,2)), cex.main = 0.75)
   }
   list(x = samples, acceptance_ratio = acceptance_ratio)
 }
 

 par(mar=c(1,1,1,1))
 unif_1 <- sim_fun(f=densFun , envelope = "unif", par1 = 0, par2 = 1, n = 10^2) 
 unif_2 <- sim_fun(f=densFun, envelope = "unif", par1 = 0, par2 = 1, n = 10^5)
 beta_1 <- sim_fun(f=densFun, envelope = "beta", par1 = 2, par2 = 2, n = 10^2)
 beta_2 <- sim_fun(f=densFun, envelope = "beta", par1 = 2, par2 = 2, n = 10^5)

 
 
 library(xtable)
 # define functions for MC estimate
 mc1 <- function(x){
   # folded normal 
   return(sqrt(2*pi/abs(x)))
 }
 
 mc2 <- function(x){
   # normal 
   exp(-x^4)/dnorm(x)
 }
 
 mc3 <- function(x){
   # uniform 
   return(20*exp(-x^4))
 }
 mcSE <- function(values){
   # this function calculates the standard error
   # of an MC estimate with target function values
   # eg values should be vector of f(x_i)
   return(sqrt(var(values)/length(values)))
 }
 # folded normal
 draws <- rnorm(100000)
 values <- sapply(draws,mc1)
 # normal
 draws2 <- rnorm(100000)
 values2 <- sapply(draws2,mc2)
 # uniform
 draws3 <- runif(100000,min=-10,max=10)
 values3 <- sapply(draws3,mc3)
 # plot histograms
 pdf(file="mc1hist.pdf")
 hist(values*2^(-5/4), xlab = expression(f(X)),
      main = "Folded Normal")
 dev.off()
 pdf(file="mc2hist.pdf")
 hist(values2, xlab = expression(f(X)),
      main = "Normal")
 dev.off()
 pdf(file="mc3hist.pdf")
 hist(values3, xlab = expression(f(X)),
      main = "Uniform")
 dev.off()
 # put results in a data frame and display with xtable
 mc.mean <- mean(values)*2^(-5/4)
 mc.sd<- 2^(-10/4)*sd(values)/sqrt(length(values))
 means <- c(gamma(1/4)/2,mc.mean,mean(values2),mean(values3))
 SEs <- c(NA, mcSE(values),mcSE(values2),mcSE(values3))
 mc.df <- data.frame(Mean = means, SE = SEs, 
                     row.names = c("True value", "row.names
                                   Normal", "Full Normal",  
                                   "Truncated Uniform"))
 xtable(mc.df)