In reikerp/contdistpeter: Peter's Continuous Distribution Package

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

This package contdistpeter currently contains two functions, rejection_sampling and d2_sampler_special.

rejection_sampling <- function(n = 1, pdf, a = 0, b = 1, C) {
  # Here I use recursion to obtain the value of a single sample.
  one_sample = function(pdf, a, b, C){
    test_draw = runif(1, a, b)
    if (runif(1, 0, C) <= pdf(test_draw)) return(test_draw)
    return(one_sample(pdf, a, b, C))
  }
  # Here I break the sampling process up using my previously defined function.
  for (i in 1:n){
    if (i == 1) my_samples = one_sample(pdf, a, b, C)
    if (i != 1) my_samples = c(my_samples, one_sample(pdf, a, b, C))
  }
  return (my_samples)
}
d2_sampler_special <- function(n = 1, fx, fyx) {
  # Determine the values to be used with the rejection_sampling function (first time).
  a = -0.01
  b = 0.01

  while(integrate(fx, -Inf, a, stop.on.error = F)$value > 1/10^4) a = a - 0.1
  while(integrate(fx, b, Inf, stop.on.error = F)$value > 1/10^4) b = b + 0.1

  # Now find an approximate value for C (first time).
  C = fx(optimize(fx, c(a,b), maximum = T)$maximum)

  # Now use rejection sampling to find a vector of sample values of x.
  my_samplesx = rejection_sampling(n, fx, a, b, C)

  # Here I initialize this vector to all 0s to save runtime.
  my_samplesy = replicate(n,0)

  for (i in 1:n) {
      # Repeat the same process to find a sample from the conditional pdf of Y given X = x.
      fyx_given_x = function(y) fyx(y,my_samplesx[i])

      j = -0.01
      k = 0.01

      while(integrate(fyx_given_x, -Inf, j, stop.on.error = F)$value > 1/10^4) j = j - 0.1
      while(integrate(fyx_given_x, k, Inf, stop.on.error = F)$value > 1/10^4) k = k + 0.1

      D = fyx_given_x(optimize(fyx_given_x, c(j,k), maximum = T)$maximum)

      y = rejection_sampling(1, fyx_given_x, j, k, D)

    my_samplesy[i] = y
  }

  my_samples = data.frame(x = my_samplesx, y = my_samplesy)

  return(my_samples)
}

Description for the rejection_sampling Function:

Single variable rejection sampling is the main function of this project. It recieves a function's pdf and generates a number of samples from the rv with the provided pdf using rejection sampling. It is thereby somewhat similar to functions such as rbinom, runif, etc.

Input variables:

n is the number of samples to generate.

pdf is the pdf of the random variable you wish to sample from.

a is the lower bound of the random variable you wish to sample from.

b is the upper bound of the random variable you wish to sample from.

C is a numeric such that f(x) <= C for all values of x.

Return variable:

my_sample A vector containg the desired number of random samples from the desired distribution.

Examples:

sample1 = rejection_sampling(10^4, function(x)1, 0, 1, 1)
sample2 = rejection_sampling(10^4, function(x)2*x, 0, 1, 2)
sample3 = rejection_sampling(10^4, function(x)exp(-x), 0, 10, 1)
hist(sample1)
hist(sample2)
hist(sample3)

Description for the d2_sampler_special Function:

This function samples from a 2D distribution given the marginal pdf of X and the conditional pdf of Y given X = x. Its output is in the form of a data frame where the elements are x and y samples from the specified 2D distribution.

This function uses the rejection_sampling function. Depending upon the sample size and complexity of the distribution, this function may run more than 1 minute.

Input variables:

n is the number of samples to generate.

fx is the marginal pdf of X.

fyx is the conditional pdf of Y given X = x (it is in the form function(y,x) expression).

Return Value:

my_samples is a data frame of x and y samples from the desired 2D distribution.

Examples:

plot(d2_sampler_special(10^3, function(x) {exp(-x^2/2)/sqrt(2*pi)}, function(y,x){exp(-y^2/2)/sqrt(2*pi)}))
plot(d2_sampler_special(10^3, function(x) {ifelse(x > 0, exp(-x), 0)}, function(y,x){ifelse(0 < y & y < x, 1/x, 0)}))
plot(d2_sampler_special(10^3, function(x) {ifelse(0.01 < x & x < 1.01, 1, 0)}, function(y,x){ifelse(0 < y, x*exp(-y*x), 0)}))