JulieLapine_Project2.R

#a) STAT 417 Project 2: Rejection Method
#b) Julie Lapine
#   19 July 2021

#c) (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.
#Create a density histogram for your generated values.
#Superimpose the probability density function of the target
#distribution to the histogram.
#Use the generated values to estimate the mean, standard deviation,
#and the probability that X < 2.)

#g) code
curve((x^2)/9, 0, 3) #graph original curve
abline(a = 1, b = 0) #uniform proposal dist - 1 is our constant

#rejection sampling below
X = runif(10000, 0, 3) #10,000 random samples from uniform dist
U = runif(10000, 0, 3)

#pi(x) function
pi_x <- function(x) {
  new_x = (x^2)/9
  return(new_x)
}

count = 1#count for loop
accept = c() #empty vector

#accept/reject algorithm
while(count <= 10000){
  test_u = U[count]
  test_x = pi_x(X[count])/(1*dunif(X[count],0,3))
  #accept conditions
  if (test_u <= test_x){
    accept = rbind(accept, X[count])
    count = count + 1
  }
  #end accept conditions
  count = count + 1
}
#end accept/reject algorithm

#d) at least part of simulated data
options(max.print = 50)
print(accept)

##               [,1]
##    [1,] 2.10430685
##    [2,] 2.98380200
##    [3,] 2.83221070
##    [4,] 2.47236733
##    [5,] 2.51535327
##    [6,] 2.93021320
##    [7,] 2.32073266
##    [8,] 2.47791505
##    [9,] 2.73889048
##   [10,] 0.88331170
##   [11,] 2.57640695
##   [12,] 2.42326209
##   [13,] 2.24272511
##   [14,] 1.97527101
##   [15,] 2.97851628
##   [16,] 2.52586890
##   [17,] 1.89460363
##   [18,] 2.80378166
##   [19,] 2.35876720
##   [20,] 2.65255862
##   [21,] 2.86111603
##   [22,] 1.94271605
##   [23,] 1.80859891
##   [24,] 2.81559059
##   [25,] 1.03495796
##   [26,] 2.09832502
##   [27,] 2.96062043
##   [28,] 2.05280095
##   [29,] 2.16371775
##   [30,] 2.36393125
##   [31,] 1.51664343
##   [32,] 1.44639493
##   [33,] 2.02817561
##   [34,] 2.17296216
##   [35,] 2.68683534
##   [36,] 1.96801829
##   [37,] 2.87507506
##   [38,] 2.97347128
##   [39,] 2.68745041
##   [40,] 1.53642013
##   [41,] 2.00835770
##   [42,] 2.93217471
##   [43,] 2.25171867
##   [44,] 2.94397241
##   [45,] 2.03751632
##   [46,] 2.22054908
##   [47,] 2.73422068
##   [48,] 2.83091093
##   [49,] 1.39784305
##   [50,] 1.32757118
##  [ reached getOption("max.print") -- omitted 2481 rows ]

#e) summary graphs/tables
hist(accept, probability = TRUE) #graph of results
curve((x^2)/9,add = TRUE, col = "red") #superimpose original pdf over results

#mean
M = mean(accept)
print(M)

## [1] 2.254159

#standard deviation
S = sd(accept)
print(S)

## [1] 0.5808469

#probability that X < 2
dnorm(2, mean = M, sd = S)

## [1] 0.6241263

#rejection rate
R = mean(1 - (accept / 10000))
print(R)

## [1] 0.9997746

#f) a discussion
#In this code, we use the rejection method to generate 10000 values from the
#equation (X^2)/9. The histogram after the fact shows how the values generally
#follow the curve of the original pdf. We also find the mean, standard deviation
#probability that x<2 from our 10,000 random variables, and print off the first
#50 values of our array.