Project 2 using the Rejection Method.

Introduction

Use the Rejection Method to generate n = 10000 values for a random variable X that has the pdf of f(x) = (x^2)/9, 0<x<3.

Part 1

Create a density histogram for your generated values.

# Declaring the Probability Density Function
f = function(x){(x^2)/9}

#Global Variables 
N = 10000
a = 0
b = 3
c = 2

count = 0 #count the number of trials 
acc.count = 0 #count the number of times Y is accepted
r = NULL # creates an empty object

for(k in 1:N)
{ 
  while(TRUE)
  {
    X = runif(1,a,b) #generates 1 random uniform value from a to b
    Y = runif(1,0,c) #generates 1 random uniform value from 0 to c
    count=count+1 #adds 1 to count every time the while loop is executed. 
    if(Y<= f(X))
    { 
      acc.count = acc.count+1 #adds 1 to acc.count if Y is less than or equal to f(X)
      break ##exits the while loop when Y is accepted
    }
  }
  r[k] = X #collect the accepted X Value 
}

head(r) #shows the first 5 values generated 

## [1] 2.8497904 2.8081825 2.8636630 2.3800847 0.7661941 2.6169833

cat("the acceptence rate is: ",acc.count/count,"\n\n")

## the acceptence rate is:  0.1652237

hist(r,prob=TRUE, xlab = 'x', main="Histogram of x")

Part 2

Superimpose the probability density function of the target distribution to the histogram.

# Declaring the Probability Density Function
f = function(x){(x^2)/9}

#Global Variables 
N = 10000
a = 0
b = 3
c = 2

count = 0 #count the number of trials 
acc.count = 0 #count the number of times Y is accepted
r = NULL # creates an empty object

for(k in 1:N)
{ 
  while(TRUE)
  {
    X = runif(1,a,b) #generates 1 random uniform value from a to b
    Y = runif(1,0,c) #generates 1 random uniform value from 0 to c
    count=count+1 #adds 1 to count every time the while loop is executed. 
    if(Y<= f(X))
    { 
      acc.count = acc.count+1 #adds 1 to acc.count if Y is less than or equal to f(X)
      break ##exits the while loop when Y is accepted
    }
  }
  r[k] = X #collect the accepted X Value 
}

#head(r) #shows the first 5 values generated 
cat("the acceptence rate is: ",acc.count/count,"\n\n")

## the acceptence rate is:  0.1660495

hist(r,prob=TRUE, xlab = 'x', main="Histogram of x")
curve(f, add=TRUE,col = "hotpink") # Theoretical probability density function

legend('topleft',
       legend = c("PDF"),
       col = c("hotpink"), 
       lwd = 2, 
       lty = c(1), 
       bty = "y")

Part 3

Use the generated values to estimate the mean, standard deviation, and the probability

that X < 2.

# Declaring the Probability Density Function
f = function(x){(x^2)/9}

#Global Variables 
N = 10000
a = 0
b = 3
c = 2

count = 0 #count the number of trials 
acc.count = 0 #count the number of times Y is accepted
r = NULL # creates an empty object

for(k in 1:N)
{ 
  while(TRUE)
  {
    X = runif(1,a,b) #generates 1 random uniform value from a to b
    Y = runif(1,0,c) #generates 1 random uniform value from 0 to c
    count=count+1 #adds 1 to count every time the while loop is executed. 
    if(Y<= f(X))
    { 
      acc.count = acc.count+1 #adds 1 to acc.count if Y is less than or equal to f(X)
      break ##exits the while loop when Y is accepted
    }
  }
  r[k] = X #collect the accepted X Value 
}

#head(r) #shows the first 5 values generated 
cat("the acceptence rate is: ",acc.count/count,"\n\n")

## the acceptence rate is:  0.1675322

hist(r,prob=TRUE, xlab = 'x', main="Histogram of x")
curve(f, add=TRUE,col = "hotpink") # Theoretical probability density function

#Simulated calculations
simulatedMean = mean(r) #calculates the average of all acceptable samples
simulatedSD = sd(r) #calculates the standard deviation of all acceptable samples 
simulatedProb = mean(r<2) #computes the proportion of values in vector r that are less than 2

#Simulated Vertical and Horizontal lines
abline(v =simulatedMean, col = "red", lty = 2, lwd = 2)
abline(v =simulatedSD, col = "orange", lty = 2, lwd = 2)
abline(h =simulatedProb, col = "yellow", lty = 2, lwd = 2)

legend('topleft',
       legend = c("SimMean","SimSD","SimProb", "PDF"),
       col = c("red","orange", "yellow","hotpink"), 
       lwd = 2, 
       lty = c(1,1,1,1), 
       bty = "y")

Part 4

Compare these estimated values with the theoretic ones.

# Declaring the Probability Density Function
f = function(x){(x^2)/9}

#Global Variables 
N = 10000
a = 0
b = 3
c = 2

count = 0 #count the number of trials 
acc.count = 0 #count the number of times Y is accepted
r = NULL # creates an empty object

for(k in 1:N)
{ 
  while(TRUE)
  {
    X = runif(1,a,b) #generates 1 random uniform value from a to b
    Y = runif(1,0,c) #generates 1 random uniform value from 0 to c
    count=count+1 #adds 1 to count every time the while loop is executed. 
    if(Y<= f(X))
    { 
      acc.count = acc.count+1 #adds 1 to acc.count if Y is less than or equal to f(X)
      break ##exits the while loop when Y is accepted
    }
  }
  r[k] = X #collect the accepted X Value 
}

#head(r) #shows the first 5 values generated 
cat("the acceptence rate is: ",acc.count/count,"\n\n")

## the acceptence rate is:  0.1703694

hist(r,prob=TRUE, xlab = 'x', main="Histogram of x")
curve(f, add=TRUE,col = "hotpink") # Theoretical probability density function

#Simulated calculations
simulatedMean = mean(r) #calculates the average of all acceptable samples
simulatedSD = sd(r) #calculates the standard deviation of all acceptable samples 
simulatedProb = mean(r<2) #computes the proportion of values in vector r that are less than 2

#Theoretical calculations
TheorMean = 2.25
TheorProb = integrate(f, lower = 0, upper = 2)$value
expectedXsquard = 5.4
TheorSD = sqrt(expectedXsquard - (TheorMean)^2)

#Simulated Vertical and Horizontal lines
abline(v =simulatedMean, col = "red", lty = 2, lwd = 2)
abline(v =simulatedSD, col = "orange", lty = 2, lwd = 2)
abline(h =simulatedProb, col = "yellow", lty = 2, lwd = 2)

#Theoretical Vertical and Horizontal lines
abline(v =TheorMean, col = "darkblue", lty = 2, lwd = 2)
abline(h =TheorProb, col = "blue", lty = 2, lwd = 2)
abline(v =TheorSD, col = "lightblue", lty = 2, lwd = 2)

legend('topleft',
       legend = c("TheorMean", "TheorProb", "TheorSD","SimMean","SimSD","SimProb", "PDF"),
       col = c("darkblue","blue", "lightblue","red","orange", "yellow","hotpink"), 
       lwd = 2, 
       lty = c(1,1,1,1,1,1,1), 
       bty = "y")

Part 5

Report the rejection rate.

# Declaring the Probability Density Function
f = function(x){(x^2)/9}

#Global Variables 
N = 10000
a = 0
b = 3
c = 2

count = 0 #count the number of trials 
acc.count = 0 #count the number of times Y is accepted
r = NULL # creates an empty object

for(k in 1:N)
{ 
  while(TRUE)
  {
    X = runif(1,a,b) #generates 1 random uniform value from a to b
    Y = runif(1,0,c) #generates 1 random uniform value from 0 to c
    count=count+1 #adds 1 to count every time the while loop is executed. 
    if(Y<= f(X))
    { 
      acc.count = acc.count+1 #adds 1 to acc.count if Y is less than or equal to f(X)
      break ##exits the while loop when Y is accepted
    }
  }
  r[k] = X #collect the accepted X Value 
}

#head(r) #shows the first 5 values generated 
cat("the acceptence rate is: ",acc.count/count,"\n\n")

## the acceptence rate is:  0.1674201

hist(r,prob=TRUE, xlab = 'x', main="Histogram of x")
curve(f, add=TRUE,col = "hotpink") # Theoretical probability density function

#Simulated calculations
simulatedMean = mean(r) #calculates the average of all acceptable samples
simulatedSD = sd(r) #calculates the standard deviation of all acceptable samples 
simulatedProb = mean(r<2) #computes the proportion of values in vector r that are less than 2

#Theoretical calculations
TheorMean = 2.25
TheorProb = integrate(f, lower = 0, upper = 2)$value
expectedXsquard = 5.4
TheorSD = sqrt(expectedXsquard - (TheorMean)^2)

#Rejection rate calculation
rejRate = 1-(acc.count/count)
cat("The Rejection Rate is:", rejRate, "\n")

## The Rejection Rate is: 0.8325799

#Simulated Vertical and Horizontal lines
abline(v =simulatedMean, col = "red", lty = 2, lwd = 2)
abline(v =simulatedSD, col = "orange", lty = 2, lwd = 2)
abline(h =simulatedProb, col = "yellow", lty = 2, lwd = 2)

#Theoretical Vertical and Horizontal lines
abline(v =TheorMean, col = "darkblue", lty = 2, lwd = 2)
abline(h =TheorProb, col = "blue", lty = 2, lwd = 2)
abline(v =TheorSD, col = "lightblue", lty = 2, lwd = 2)

legend('topleft',
       legend = c("TheorMean", "TheorProb", "TheorSD","SimMean","SimSD","SimProb", "PDF"),
       col = c("darkblue","blue", "lightblue","red","orange", "yellow","hotpink"), 
       lwd = 2, 
       lty = c(1,1,1,1,1,1,1), 
       bty = "y")

Discussion

The rejection method allowed for the simulation of a PDF function that is not uniform. the rejection method created a way to generate samples that followed the f(x) distribution by using the uniform random number. The method accepted values based on the probability that they would be under the desired PDF.