knitr::opts_chunk$set(echo = TRUE)
suppressPackageStartupMessages(library(tidyverse))

Part a:

given that K = sup(f(x)/g(x)), show that K = sqrt(2e/pi)

function f(x)/g(x) plotted below

fx <- function(x) {
  (1/sqrt(2*pi))*exp((-x^2)/2)
  }
gx <- function(x) {
  (1/2)*exp(-abs(x))
}


pdfs <- data.frame(x = runif(n = 15000, min = -10, max = 10)) %>%  
  mutate(fx = fx(x), 
         gx = gx(x), 
         fxgx = fx/gx)

ggplot(pdfs)+
  geom_point(aes(x, fxgx), color = "blue")+
  labs(y = "f(x)/g(x)")

the x-values that the suprema occur at are -1 and 1, you could also show this by taking the derivative of f(x)/g(x) and setting it equal to zero to find local maxima.

pdfs$x[pdfs$fxgx == max(pdfs$fxgx)]
## [1] -0.9999803

we can then set our equations f(x) and g(x) to have these values

max1 <- fx(1)/gx(1)
max2 <- fx(-1)/gx(-1)


maxes <- c(max1, max2)

k = sqrt((2*exp(1))/pi)

all(maxes == k) #all proofs equal to K? yes
## [1] TRUE

Part B plot f(x) and Kg(x) on the same plot

pdfs <- pdfs %>%  
  mutate(kgx = gx*k)

ggplot(data = pdfs)+
  geom_point(aes(x, fx), color = "blue")+ 
  geom_point(aes(x, kgx), color = "red")+
  labs(caption = "red line is K * g(x)\nblue line is f(x)")+ 
  labs(y = "Output")

## c, d: try to do this on your own so you can understand what’s going on with this

part e implement the algorithm

accept_reject <- function(n) {
  
  pdf <- data.frame(x = runif(n = n, min = -10, max = 10), y = runif(min = 0, max = 0.5, n = n)) %>%  
    mutate(gx = gx(x), 
           accept = ifelse(y < gx, y, NA))
  accepts <- sum(!is.na(pdf$accept))
  rejects <- sum(is.na(pdf$accept))

  pdf <- pdf %>% filter(!is.na(accept)) #removes all of the rejects
  
  while(nrow(pdf) < n-(0.1*n)) {#semi-vectorized function up until about 90% of accepts are met
    
    temp <- data.frame(x = runif(n = n, min = -10, max = 10), y = runif(min = 0, max = 0.5, n = n)) %>%  
    mutate(gx = gx(x), 
           accept = ifelse(y < gx, y, NA))
    
    accepts <- accepts + sum(!is.na(temp$accept))
    rejects <- rejects + sum(is.na(temp$accept))
    
    temp <- temp %>% filter(!is.na(accept))
    pdf <- rbind(pdf, temp) 
    
  }
  
  while(nrow(pdf) < n) { #one-by-one while loop for the last ~10% of accepts
    temp <- data.frame(x = runif(n = 1, min = -10, max = 10), y = runif(min = 0, max = 0.5, n = 1)) %>%  
    mutate(gx = gx(x), 
           accept = ifelse(y < gx, y, NA))
    
    accepts <- accepts + sum(!is.na(temp$accept))
    rejects <- rejects + sum(is.na(temp$accept))
    
    temp <- temp %>% filter(!is.na(accept))
    pdf <- rbind(pdf, temp)
    
  }
  
  efficiency <- paste0("Efficiency: ", accepts/(accepts+rejects))
  print.noquote(efficiency)
  
  pdf <- pdf %>% 
    select(x, accept)
  return(pdf)
}


final <- accept_reject(5000)
## [1] Efficiency: 0.0982356870603953
ggplot()+
  geom_point(data = final, aes(x, accept), alpha = 0.25)+
  geom_point(data = pdfs, aes(x, gx), color = "red", alpha = 0.1)+
  geom_point(data = pdfs, aes(x, fx), color = "blue", alpha = 0.1)+
  labs(caption = "Red is the Double Standard Equation\nBlack are the points obtained from the accept/reject algorithm using double standard equation as proposed pdf\nBlue is the normal distribution", title = "Accept/Reject Algorithm Using Double Standard Equation")

ggplot()+
  geom_histogram(data = final, aes(x))+
  labs(title = "Accepted Values Using Accept/Reject Algorithm", 
       caption = "n=5000")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

How well does this proposed pdf fit to the normal distribution?

Mean of this simulated distriution using N = 5000 is -0.0451041. Standard deviation of this distribution is 1.4026997.