In bubj/samplr:

knitr::opts_chunk$set(
  collapse = TRUE
)
library(samplr)
library(ggplot2)

This package implements one and two dimensional rejection sampling for continuous probability distribution functions.

One-dimensional

The easiest form of rejection sampling occurs using one-dimensional pdfs. For example, we can sample from a beta distribution, in this case Beta(2,1). We first define a function for this pdf, where the varaible must be defined as $x$.

betapdf <- function(x) {
  ifelse(0 < x & x < 1, 2*x, 0)
}

We can then use the samplr function to sample from this distribution,

samplr(1,betapdf)

With many samples, we are able estimate the distribution by using hist to plot the density of the samples,

hist(samplr(10000,betapdf), main = NULL)

We can also sample from distributions with unbounded support,

f <- function(x) {
  if (x < 0) {
    1/2*exp(x)
  }
  else {
    1/2*exp(-x)
  }
}
hist(samplr(10000,f), main = NULL)

We can find probabilities associated with these pdfs using the psamplr function, which is similar to the built in R functions pnorm, pbinom, etc., which gives $P(X < q)$, for some $q$ specified in the calling of the psamplr function.

psamplr(betapdf,.5)
psamplr(f,0)

Expectation can also be calculated using the expsamplr function, which calculates $E[X]$.

expsamplr(betapdf)
expsamplr(f)

Two-dimensional

We can also use the samplr function to sample from two-dimensional pdfs as well, specifying the parameter twod = TRUE. For two-dimensional pdfs, the parameter of the function that you wish to sample from must be a vector containing the random variables $x$ and $y$. The output of the function is a matrix of N rows, where each row corresponds to a pair of points, $x$ and $y$.

g <- function(z) {
  x <- z[1]
  y <- z[2]
  ifelse(0 < x & x < 1 & 0 < y & y < 1, x + y, 0)
}
samplr(1, g, twod = TRUE)

We can get an idea of what this distribution looks like by finding many samples, and using ggplot with geom_density_2d() to get contours of the pdf.

samps <- data.frame(samplr(10000, g, twod = TRUE))
colnames(samps) <- c("x","y")
ggplot(samps, aes(x, y)) + geom_density_2d()

The samplr function also works for two-dimensional pdfs with unbounded support,

h <- function(z) {
  x <- z[1]
  y <- z[2]
  ifelse(0 < x & 0 < y, exp(-x)*exp(-y), 0)
}
samps <- data.frame(samplr(10000, h, twod = TRUE))
colnames(samps) <- c("x","y")
ggplot(samps, aes(x, y)) + geom_density_2d()

Like for one-dimensional pdfs, probabilities can also be calculated for two-dimensional pdfs using the psamplr function with a second function used to take various combinations of the random variables X and Y. The second function, g in the parameters, must have an input that is a vector of the random variables. For two-dimensional pdfs, psamplr calculates $P(g(X,Y) < q)$.

prob1 <- function(z) {
  x <- z[1]
  y <- z[2]
  x + y
}
prob2 <- function(z) {
  x <- z[1]
  y <- z[2]
  x*y
}
psamplr(g,1,prob1)
psamplr(g,1,prob2)

Expectation values can also be found for two-dimensional pdfs using the expsamplr with a second input function of a combination of the random varialbes, similar to the psamplr function. For two-dimensional pdfs, expsamplr calculates $E[g(X,Y)]$.

expectx <- function(z) {
  x <- z[1]
  y <- z[2]
  x
}
expsamplr(g,prob1)
expsamplr(g,expectx)

bubj/samplr documentation built on May 28, 2019, 7:14 p.m.