The Problem

What to Expect

The rejection method is a sampling method to deal with hard to sample functions. To deal with the hard function, we will generate a random numbers from a uniform distribution and accept those that fall within the function and reject those that do not.

Since we are looking at the function from \(x=0\) to \(x=3\), we must find what the min and max values are for the function in that range.

The function we are looking at is an exponential function, and in this case, we will not need to worry about maximums or minimums that are not obvious.

f=function(x) (x^2) / 9
x.lower = 0
x.upper = 3
f.min = f(x.lower)
f.max = f(x.upper)

Here, we find that our minimum \(f(0)=0\) and our maximum \(f(3)=1\).

Since the maximum is 1, that is what we will set our \(c\) value to. Our \(a\) and \(b\) values will match our minimum and maximum \(x\) values, respectively.

c = f.max
a = x.lower
b = x.upper

We will calculate the expected mean, standard deviation, and probability of \(X<2\) for comparison later.

# Find the expected mean
f.mu = integrate(function (x) {x * f(x)},
                 lower = x.lower,
                 upper = x.upper)$value
# Find the expected standard deviation
f.sd = (integrate(function (x) {x^2 * f(x)},
                  lower = x.lower,
                  upper = x.upper)$value - f.mu^2)^(1/2)
# Find the probability X<2
f.p2 = integrate(f,
                 lower = x.lower,
                 upper = 2)$value

The Work

With our parameters set, we can setup the simulation and start to determine our accept and reject numbers.

# Number of simulations
n = 10000
accept = NULL
accept.count = 0
i = 0

for(i in i:n){
  X = runif(1,
            min=a,
            max=b)
  Y = runif(1,
            min=0,
            max=c)
  if (Y <= f(X)){
    accept.count = accept.count+1
    accept[accept.count] = X
  } 
}

The following displays the first 50 values of the simulated data.

head(accept, n = 50)
##  [1] 2.980112 2.205976 1.282570 1.747612 2.135702 2.709253 1.864181 2.861321
##  [9] 1.993022 2.904558 1.441091 1.436229 2.904970 1.402810 1.922313 2.739116
## [17] 2.593193 1.570468 2.982636 2.316899 2.553509 2.642389 2.457233 1.951445
## [25] 1.753927 1.472277 1.427664 2.376850 2.896913 1.628960 2.400527 2.946566
## [33] 2.262862 2.738661 1.181250 2.827985 2.304734 1.648482 2.947182 1.776121
## [41] 2.493352 2.249762 2.683493 1.271735 2.882608 1.883727 2.542684 1.187195
## [49] 1.951839 2.717928

Presenting the Data

hist(accept,
     probability =  TRUE)
curve(f(x),
      add = TRUE)

Findings

calc.mu = mean(accept)
calc.sd = sd(accept)
# To find the probability X<2, count the number accepted < 2 / all accepted
calc.p2 = sum(accept < 2) / accept.count

# Rejection rate
1 - (accept.count / n)
## [1] 0.6609

The simulated data matches up closely with our expected data. We expected a mean of 2.25 and received a calculated mean of 2.217, with a difference of 0.0327.

Our standard deviation we expected was 0.5809 and our data showed a standard deviation of 0.598, with a difference of 0.017.

Additionally, we calculated the probability of \(X<2\) as 0.2963, whereas our count gave the probability as 0.3144, with a difference of 0.0181

Conclusions

From this simulation, we find that, the rejection method, can be very useful to run a simulation on a function that may not be easy to simulate on its own. By using a random uniform distribution and only accepting those numbers that match our criteria, we calculate statistics for some very difficult formulae.