simulation-1

Overview of Topics on Simulation

• Exact Simulation

  • Simulating from Standard Parametric Distributions

  • Quantile transform method (Inverse CDF)

• Monte Carlo Simulation

  • Monte Carlo integration

• Markov Chain Monte Carlo (MCMC)

  • Metropolis-Hastings Algorithm

  • Gibbs Sampling

  • Validate Computation: Convergence diagnostics, Simulation-based calibration (optional)

  • Output analysis for MCMC

Introduction

Why simulate?

• We want to see what a probability model actually does
• We want to understand how our procedure works on a test case
• We want to use a partly-random procedure

All of these require drawing random variables from distributions

We have seen R has built in distributions: beta, binom, cauchy, chisq, exp, f, gamma, geom, nbinom, norm, pois, t, unif, etc.

Every distribution that R handles has four functions:

• p for “probability”, the cumulative distribution function (c. d. f.)

• q for “quantile”, the inverse c. d. f.

• d for “density”, the density function (p. f. or p. d. f.)

• r for “random”, a random variable having the specified distribution

#   Densities for various distributions
set.seed(1234)
dnorm(0, mean=0, sd=1)      # Normal distribution
dexp(1, rate=1)         # Exponential distribution
dbinom(5, size=10, prob=0.5)    # Binomial distribution
dpois(10, lambda=10)        # Poisson distribution
dgamma(1, shape=1, scale=1) # Gamma distribution

## NOTE: Here the density function is being evaluated at the given x with the specified parameters

# CDFs, inverse CDFs, and random generation

pnorm(0, mean=0, sd=1)      # P(X<=0)
qnorm(0.95, mean=0, sd=1)   # quantile related to 95% lower tail probability
rnorm(10, mean=0, sd=1)     # generates 10 realizations of a standard normal RV

Exact Simulation

Simulating from Standard Parametric Distributions

• Usually, R gets Uniform\((0,1)\) random variates via a pseudo-random generator
• R uses a sequence of Uniform\((0,1)\) random variates to generate other distributions

Example: Binomial

• Suppose we want to generate a Binomial\((1,1/3)\) using a \(U \sim \text{Uniform}(0,1)\)
• Consider the function \(X^* = I(0<u<1/3)\), then \[ P(X^* = 1) = P( I(0<u<1/3) = 1) = P( u \in (0,1/3)) = 1/3 \] and \(P(X^* = 0) = 2/3\)
• Hence, \(X^* \sim \text{Binomial}(1,1/3)\)
• Two ways to extend this to Binomial\((n,1/3)\)

my.binom.1 <- function(n=1, p=1/3){
  u <- runif(n)
  binom <- sum(u<p)
  return(binom)
}

my.binom.1(1000)

[1] 351

my.binom.1(1000, .5)

[1] 465

my.binom.2 <- function(n=1, p=1/3){
  u <- runif(1)
  binom <- qbinom(u, size=n, prob=p)
  return(binom)
}

my.binom.2(1000)

[1] 341

my.binom.2(1000, .5)

[1] 515

Quantile transform method (Inverse CDF)

Given \(U \sim \text{Uniform}(0,1)\) and CDF \(F\) from a continuous distribution. Then \(X = F^{-1}(U)\) is a random variable with CDF \(F\). \[ P(X\le a) = P (F^{-1}(U) \le a) = P ( U \le F(a)) = F(a) \]

• \(F^{-1}\) is the quantile function
• If we can generate uniforms and calculate quantiles, we can generate non-uniforms
• Also known as the Probability Integral Transform Method

Example: Exponential

Suppose \(X \sim \text{Exp}(\beta)\). Then we have density \[ f(x) = \beta^{-1} e^{-x/\beta} I(0<x<\infty) \] and CDF \[ F(x) = 1 - e^{-x/\beta} \] Also \[ u = 1 - e^{-x/\beta} \text{ iff } -x/\beta = \log (1-u) \text{ iff } x = -\beta \log (1-u). \] Thus, \(F^{-1} (u) = -\beta \log(1-u)\).

So if \(U \sim \text{Uniform}(0,1)\), then \(F^{-1} (u) = -\beta \log(1-u) \sim \text{Exp}(\beta)\).

x <- runif(10000)
y <- - 3 * log(1-x)
hist(y)

mean(y)

[1] 3.013837

true.x <- seq(0,30, .5)
true.y <- dexp(true.x, 1/3)
hist(y, freq=F, breaks=30)
points(true.x, true.y, type="l", col="red", lw=2)