Using the Rejection Method

  1. Use the Rejection Method to generate \(n = 10000\) values for a random variable X that has the pdf of \[f(x) = \frac{x^2}{9}~~,~~0<x<3\]
  2. Create a density histogram for your generated values.
  3. Superimpose the probability density function of the target distribution to the histogram.
  4. Use the generated values to estimate the mean, standard deviation, and the probability that \(X<2\).

Rejection Method

We will use the uniform distribution which has the same support as the desired distribution: \(Uniform(0,1)\).

We begin by creating our function and generating \(n\) values:

N = 10000     # number fo values
i = 0         # number of candidate values before N valid values are generated using the proposal distribution
n = 0         # iterations
y = NULL      # empty vector
# our function
f <- function(x) {
  fx = (x**2)/9
  return(fx)
}

# generating values
while (n <= N) {
  i = i + 1             # iterator
  x = runif(1, 0, 3)    # random number generated
  u = runif(1)          # uniformly generated value
  c = 3                 # max value of ratio
  a = (9*f(x))*(1/c)    # alpha value used for rejection
  
  # accepting/rejecting values
  if (u <= a) {
    n = n + 1
    y[n] = x
  }
}

# A little sneak at our values
y[1:50]
##  [1] 2.8004782 2.7856484 1.9038437 2.2662282 2.6500598 1.4019888 2.9305635
##  [8] 2.3852936 2.3571833 2.1281242 2.4932133 1.7717072 0.9537013 2.1117129
## [15] 2.0275985 2.0351414 1.5040497 2.5533518 1.3864536 1.8112198 2.9859133
## [22] 1.4610562 0.7380692 1.1967870 2.4992792 2.2290436 0.7121019 1.7471602
## [29] 2.6671671 0.4975265 1.9645273 2.8337730 1.6766634 1.2707174 2.4294154
## [36] 1.2403571 1.2878072 2.7734655 1.2704538 2.4042130 1.9645495 1.6441994
## [43] 1.1977579 2.7821267 2.1054670 2.4338237 2.1922303 1.9642233 2.5718712
## [50] 1.8039446

Now we create a visual for our generated distribution and superimpose a probability density function of our target distribution.

hist(y, probability = TRUE, xlab = "x", main = "Histogram of the Generated Sample")
curve(f, add = TRUE, lty = "dashed", col = "red")

cat("Our rejection rate is",1-(N/i),".")
## Our rejection rate is 0.3857117 .

Using Our Distribution

m = mean(y)
stddev = sd(y)

cat("Our distribution has the mean", m, "and the standard deviation", stddev)
## Our distribution has the mean 2.046503 and the standard deviation 0.6086111

The probability that \(X<2\) is:

pexp(2, 1/m)
## [1] 0.6236655

Summary

Our function was expected to have an underlying exponential distribution since it is a power function. We used the rejection method to generate our own distribution for that specific function (our target distribution).