Original Description

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.

  1. Create a density histogram for your generated values.
  2. Superimpose the probability density function of the target distribution to the histogram.
  3. Use the generated values to estimate the mean, standard deviation, and the probability that X < 2.
  4. Compare these estimated values with the theoretic ones.
  5. Report the rejection rate.

Code

f=function(x) {
  return((x^2)/9)
} #function is the pdf of the target distribution
N=10000 #size of sample to generate
y=numeric(N) #initial states of N values
i=0 #holds values
n=0 #controls iterations
while(n<=N){
  i=i+1
  x=runif(1,0,3) #distribution is uniform (0,3)
  u=runif(1)
  c=1 #maximum value of ratio (f(x)=g(x))
  alpha=f(x)/((3^2)/9)/c
  if (u<=alpha){
    n=n+1 #accepted values increases by 1
    y[n]=x #stores value
  }
}

Part of Simulated Data

y[1:20]
##  [1] 2.746691 2.876025 2.847241 2.299315 1.903207 2.455342 2.607857 2.794663
##  [9] 2.244020 2.463646 2.339030 1.646711 1.695895 2.603744 2.753719 1.017442
## [17] 2.086237 2.934437 2.972996 1.262959

Summary Graph

hist(y, probability=TRUE, xlab="x", main="Histogram")
#superimpose pdf on the histogram
curve(f, add=TRUE, lty="dashed", col="red")
abline(v = mean(y),   col = 2, lty = 2, lwd = 2)
legend(2, 0.5, c("Simulated Mean"), col = 2:3,  text.col = 2:3, lty = 2:3, lwd = 2, bty = "N")

The Rejection Rate, Mean, Standard Deviation, and Probability that X<2

cat("The rejection rate: ", 1-(N/i), ".")
## The rejection rate:  0.6664443 .
cat("The acceptance rate: ", N/i, ".")
## The acceptance rate:  0.3335557 .
cat("The simulated mean: ", mean(y), ".")
## The simulated mean:  2.259585 .
cat("The theoretical mean: ", 2.25, "." )   #prints the theoretical mean
## The theoretical mean:  2.25 .
cat("The simulated probability that X<2: ", mean(y<2), ".") #estimates the simulated probability X<2
## The simulated probability that X<2:  0.2918708 .
cat("The theoretical probability that X<2: ", 0.296, ".") #calculates the theoretical probability X<2
## The theoretical probability that X<2:  0.296 .
cat("The simulated standard deviation: ", sd(y), ".") #estimates the simulated std dev
## The simulated standard deviation:  0.5798414 .
cat("The theoretical standard deviation: ", 0.58094750, ".")#prints the theoretical std dev
## The theoretical standard deviation:  0.5809475 .

Discussion

In this project, we generated 10000 values and a density histogram with the probability density function superimposed using the rejection method. Additionally, we generated the rejection rate, acceptance rate, mean, standard deviation, and probability that X<2. The rejection rate was calculated to be approximately 0.6658 and the acceptance rate was found to be around 0.3342. The other simulated results can be compared to their theoretical values that were calculated separately from the code, respectively. First, the theoretical mean was 2.25, while our simulated mean was approximately 2.2482, which are quite close in value. Then, the theoretical standard deviation was found to be around 0.5809 with the generated standard deviation resulting in around 0.5858. The two values compare to be somewhat close in value, but not as equivalent with approximation as the mean. Finally, the theoretical and simulated probability that X<2 showed to be similar to the values associated to the mean, that is, being almost equivalent, with resulting approximate values being 0.296 and 0.2967, respectively.