Simulation

Simulation Methods

• From a uniform random variable \(U\sim Uni(0,1)\) we can generate a RV from any arbitrary distribution or approximate integrals. The uniform random variable, with density defined as follows, is our source of randomness for further simulations.

\[ f(u) = \begin{cases}1, \quad\text{if }u\in[0,1]\\0,\quad\text{otherwise}\end{cases}\]

Linear Congruential Method For Uniform Simulation

• we can simulate a sequence of random numbers uniform on \([0,1]\) through the following method:

  • Define \(0\leq x_0\leq m-1\), as well as constants \(a,b,m\). These constants need to be chosen carefully so that \(u_i\) shows characteristics of independent random uniform variables.

  • Define sequence \(\lbrace x_n:n\geq 0\rbrace\) by recursion, noting that

    \[ x_{k+1}=(ax_k+b)\text{ mod } m\]

  • We approximate sequence of random uniform numbers by the following formula:

    \[ \lbrace u_n = \frac{x_n}{m}\rbrace,\quad 0\leq u\leq 1 \]

Approximating Integrals

• Using the random uniform numbers we can approximate integrals in the following form numerically:

\[ \int^1_0g(u) du\]

• As \(U\sim Uni(0,1)\), \(f(u)=1\in[0,1]\), therefore we can use the following approximation:
  • As \(n\rightarrow\infty\) the approximation converges to the real value

\[ \int^1_0g(u)f(u) du=\int^1_0g(u)du = \mathbb{E}[g(U)]\]

\[ \mathbb{E}[g(U)]\approx\frac{\sum^n_{i=1}g(U)}{n} \]

• In the case where \(a,b\) are finite, we scale to \(0,1\) by using the transformation \(u=\frac{(x-a)}{(b-a)}\), then \(x=a+u(b-a)\):

\[ \int^b_a g(x) dx = \int^1_0(b-a)g(a+u(b-a)) du\]

• Suppose range is infinite, we make the transformation \(u=\frac{1}{(x+1)}\) such that \(x=\frac{1}{u}-1\):
  • If one of the bounds are a constant \(c\), we use \(\frac{1}{x+1-(c)}\)
  • If the lower bound is \(-\infty\) use \(u = \frac{1}{1-x}\) such that \(x = 1-\frac{1}{u}\)

\[ \int^\infty_0 g(x) dx = \int^1_0\frac{1}{u^2}g\left(\frac{1}{u}-1\right) du\]

• Monte Carlo methods are useful for approximation of multidimensional integrals aswell
  • This can be generalised for higher dimensions, but can become inefficient as dimensions increase as the sam-ples required increases dramatically

\[ \int^1_0\int^1_0g(u_1,u_2)f(u_1,u_2) du_1du_2 = \int^1_0\int^1_0g(u_1,u_2) = \mathbb{E}[g(U_1,U_2)]\approx\frac{\sum^n_{i=1}g(U_1^i,U_2^i)}{n}\]

• A simple example in R is as following:

\[ \int^2_0\exp(x^2)dx = 2\int^1_0\exp(4u^2) du\]

u = runif(1000)
mean(2*exp(4*u^2))

## [1] 16.36616

• A more complicated example:

\[ \int^\infty_{-2}\exp(-x^2) dx= \int^1_0\frac{1}{u^2}\exp\left(-\left(\frac{1}{u}-3\right)^2\right)du\]

u=runif(100000)
mean((1/(u^2))*exp(-((1/u)-3)^2))

## [1] 1.771085

Generating Discrete Random Variables

Inversion Method

• We can consider \(Y\) as a random variable with probability function \(\lbrace p_j\rbrace\):

\[ P(Y=j)=P\left(\sum^{j-1}_{i=0}p_i\leq U < \sum^j_{i=0}p_i\right) = \left(\sum^j_{i=0}p_i\right)-\left(\sum^{j-1}_{i=0}p_i\right)= p_j\]

EDF

• An example in R:

\[p_j = P(Y=J)=\exp(-\lambda)\frac{\lambda^j}{j!}\]

• Note that in this case we can use the recursive relationship to make things easier:

\[ p_{j+1} = \exp(-\lambda)\frac{\lambda^{j+1}}{j+1!}\quad\therefore\quad p_{i+1}=\frac{\lambda}{i+1}p_i\]

• In the following R code, we are checking each step and seeing if the uniform random variable falls within each

lambda = 2
u = runif(1)
p <- exp(-lambda) # Setting p_0
cumsum <- p
y <- 0
while (u>cumsum){
  p<- lambda/(y+1)*p
  cumsum <- cumsum+p
  y=y+1
}

Rejection Sampling

• The inversion technique isn’t very general and can be hard or impossible to implement for some distributions. Therefore we can also consider the rejection method, which is more general and therefore works for more distributions.
• We consider \(p_j,\quad j=0,1,...\), the distribution we want to simulate from and \(q_j,\quad j=0,1,...\) an easy to simulate from distribution.
• We assume there exists a constant \(c\) such that:
  • This \(c\) is the maximum ratio

\[ \frac{p_j}{q_j}\leq c\]

• The rejection algorithm is as follows:
  • Simulate a random variable \(Y\) from probability mass function \(\lbrace q_j\rbrace\)
  • Generate a uniform random variable \(U\) on \([0,1]\)
  • If \(U<\frac{p_Y}{cq_Y}\) then let \(X=Y\) and stop. Otherwise, go to step 1.
• This works as

\[ P(X=j) = P(Y=j|Accept) = \frac{p_j}{cK},\quad\therefore\quad cK=1\text{ and }P(X=j)=p_j\]

• From this we can see that the acceptance probability is \(K=\frac{1}{c}\)
  • Therefore as \(c\rightarrow\infty\), \(K\rightarrow 0\) and the rejection sampling method becomes increasingly inefficient.
• A simple example in R
  • We want to simulate a RV \(X\) taking one value \(1,...,10\) with respective probabilities stated in the code. We use the rejection method using a uniform distribution on the integers \(1,...,10\), aka \((q_j=0.1)\)

p <- c(0.11,0.12,0.09,0.08,0.12,
       0.1,0.09,0.09,0.1,0.10)
q <- rep(0.1,10)
c <- max(p/q)
c

## [1] 1.2

y <- sample(1:10,1)
while(runif(1)>p[y]/(0.1*c)){
  y <- sample(1:10,1)
}
y

## [1] 1

Generating Continuous Random Variables

Inverse Transform Method

• Simulation of a continuous random variable with the given distribution function \(F(x)\) can be done with the use of a uniformly distributed random variable. We let \(F(x)\) be the distribution function, supposing it is invertible, we write \(F^{-1}(p)\) for the inverse. Therefore, we can simulate \(X = F^{-1}(U)\)

• If \(F(x)\) is not invertible there is a general result which states that if we define \(F^-(p)=inf\lbrace x:F(x)\geq p\rbrace\) then \(F^-(U)\) has distribution function \(F(x)\).

• Note that for an inverse:

\[ F(F^{-1}(p))=p\]

• Example with R Code:
• Suppose \(Y\sim \exp(\lambda)\)

\[ f(y) = \begin{cases}\frac{1}{\lambda}\exp(-\frac{y}{\lambda}),\quad y\geq 0\\ 0,\quad\text{otherwise}\end{cases} \]

\[ F(y) = \begin{cases}1-\exp(-\frac{y}{\lambda}),\quad y\geq 0\\ 0,\quad\text{otherwise}\end{cases}\]

\[ 1-\exp(-F^{-1}(p)/\lambda)=p,\quad \therefore\quad F^{-1}(p) = -\lambda\log(1-p)\] * A R function that simulates this variable would therefore be: + Noting \(\log(U)\) is the same as \(\log(1-U)\) as \(U\) and \(1-U\) are both uniform on \([0,1]\)

simexp <- function(n,lambda){
  u <- runif(n)
  x <- (-lambda)*log(u)
  x
}
simexp(10,3)

##  [1]  0.8453050  0.8162995  1.2715004  2.6992083  3.7324242  2.2712537
##  [7]  1.3701184  1.7007608  1.2672142 10.9796646

• A more complicated example:

\[ f(x) = \begin{cases}\frac{(x-2)}{2},\quad \text{if }2\leq x < 3\\ \frac{(2-x/3)}{2},\quad\text{if }3\leq x\leq 6\\ 0,\quad\text{otherwise}\end{cases}\]

\[ x<2, F(x)=0\]

\[ 2\leq x<3,\quad F(x)=\int^x_2\left(\frac{z-2}{2}\right)dz=\frac{(x-2)^2}{4}\]

\[ 3\leq x\leq 6,\quad F(x)=\int^x_3\left(\frac{2-z/3}{2}\right)dz=1-\frac{3}{4}(2-x/3)^2\] \[ x>6, F(x) = 1\]

• Therefore When \(F^{-1}(p)\in(2,3]\), aka \(p\in(F(2),F(3)]=(0,1/4])\)

\[ (F^{-1}(p)-2)^2/4=p,\quad\therefore F^{-1}(p)=2+\sqrt{4p} \]

• Therefore When \(F^{-1}(p)\in(3,6]\), aka \(p\in(F(3),F(6)]=(1/4,1])\)

\[ 1-3/4(2-F^{-1}(p))^2=p,\quad\therefore F^{-1}(p)=3(2-\sqrt{4/3(1-p)})\]

simfx <- function(n){
  u <- runif(n)
  x <- ifelse(u<=1/4, 2+sqrt(4*u), 3*(2-sqrt(4/3*(1-u))))
  x
}
simfx(10)

##  [1] 4.066587 4.002064 3.579809 2.165836 4.847225 3.079818 3.813328 3.369304
##  [9] 4.357349 2.673095

Rejection Sampling

• Similar idea to the discrete case however we used density functions:

\[ \frac{f(x)}{g(x)}\leq c\]

• Simulate with the algorithm:
  • Generate \(Y\) with density \(g\)
  • Generate \(U\sim Uni(0,1)\)
  • If \(U\leq\frac{f(Y)}{cg(Y)}\) set \(X=Y\). Otherwise go to step 1.
• Example with R code:

\[ f(x) =\begin{cases}20x(1-x)^3,\quad, 0\leq x\leq1\\ 0,\quad\text{otherwise}\end{cases}\]

\[ g(x) = Uni(0,1) \]

• Find \(c\) by maximising the ratio, which in the case of uniform is maximising \(f(x)\)
  • Global maximum (check second derivative) occurs at \(x=1/4\), aka \(f(1/4)=135/64=c\)

\[ f'(x) = 20(1-x)^2(1-4x)\]

\[ \frac{f(x)}{g(x)}\leq\frac{135}{64}\]

simreject <- function(){
  y <- runif(1)
  c <- 135/64
  while (runif(1)>(20*y*(1-y)^3)/c){
    y<- runif(1)
  }
  x<- y
  x
}
simreject()

## [1] 0.1885538

Analysis Of Simulated Data

• We often use simulation to estimate quantities of interest for some stochastic system.
• Estimate of \(\theta\): The sample mean

\[ \bar{X}=\frac{1}{n}\sum^n_{i=1}X_i,\quad \mathbb{E}[\bar{X}]=\theta,\quad\mathbb{V}ar(\bar{X}) = \frac{\sigma^2}{n}\]

• We estimate \(\sigma^2\) by taking the sample variance:

\[ S^2 = \frac{\sum^n_{i=1}(X-\bar{X})^2}{n-1}\]

• \(\frac{S^2}{n}\) is therefore our estimate for the variance.

• We are often interested in how many samples, \(n\), we should generate for a specific task. We often implement a stopping procedure to choose in this case:

  • Choose an acceptable value \(d\) for standard deviation of the estimator
  • Generate at least 30 data values
  • Continue to generate data values, stopping when you have generated \(k\) values and \(S/\sqrt{k}\)<d

• For efficiency, we create a recursive formula to calculate the mean and variance.

\[ \bar{X}_{j+1}=\bar{X}_j+\frac{X_{j+1}-\bar{X}_j}{j+1} \]

\[ S^2_{j+1}=\left(1-\frac{1}{j}\right)S^2_j+(j+1)(\bar{X}_{j+1}-\bar{X}_j)^2\]

• Example with R Code:
  • We want to estimate the following value when the estimated standard deviation \(d\) is less than 0.01

\[ \theta = \int^1_0\exp(x^2)dx\]

u <- runif(100)
cmean <- mean(exp(u^2))
n <- 100
d <- 0.01
s <- sqrt(var(exp(u^2)))
while(s/sqrt(n)>d){
  newx <- exp(runif(1)^2)
  nmean <- cmean+(newx-cmean)/(n+1)
  s <- sqrt((1-(1/n))*s^2+(n+1)*(nmean-cmean)^2)
  cmean <- nmean
  n <- n+1
}
cmean

## [1] 1.478606

n

## [1] 2329

• We may also want to give an interval estimate of \(\theta\).
  • We can use the central limit theorem to get such an estimate

\[ \frac{\bar{X}-\theta}{\sigma/\sqrt{n}}\rightarrow N(0,1)\]

• Therefore an approximate \(100(1-\alpha)\) percentage confidence interval for \(\theta\) is:

\[ \bar{X}\pm z_{\alpha/2}\frac{\sigma}{\sqrt{n}}\]

• We can base our decision on when to stop simulation on achieving a certain width for a confidence interval rather than a given sd.
  • The width of a \(100(1-\alpha)\%\) interval for \(\theta\) is \(2z_{\alpha/2}s/\sqrt{n}\)
  • Choose an acceptable \(l\) for the width of the confidence interval
  • Generate atleast 30 data values
  • Continue generation until \(2z_{\alpha/2}S/\sqrt{k}<l\)
• Consider the same problem as before, however we want a 95% confidence interval to have width \(<0.01\)

u <- runif(100)
cmean <- mean(exp(u^2))
n <- 100
s <- sqrt(var(exp(u^2)))
while (2*1.96*s/sqrt(n)>0.01){
  newx <- exp(runif(1)^2)
  nmean <- cmean + (newx-cmean)/(n+1)
  s <- sqrt((1-1/n)*s^2+(n+1)*(nmean-cmean)^2)
  cmean <- nmean
  n <- n+1
}
cmean

## [1] 1.457839

n

## [1] 34067

Bootstrap

• Bootstrap is useful when we have an estimator of parameter \(\theta\) that does not take the form of a sample mean. The bootstrap approach is a very general method for assessing the properties of an estimator \(\hat{\theta}\) of \(\theta\)

• The idea of bootstrap is to sample the empirical distribution function in order to estimate quantities of an unknown underlying \(F(x)\).

• For bootstrap we sample a value uniformly from the original sample with replacement, which is the same as simulating from the empirical distribution function \(\hat{F}\)

• Example using the variance of the sample mean to show how it works (bootstrap is not needed for sample mean variance)

• Bootstrap approximation of sample mean:

\[ \mathbb{E}(Y)=\sum^n_{i=1}x_i\underbrace{P(Y=x_i)}_{1/n}=\bar{x}\]

• Bootstrap approximation of sample variance:

\[ \mathbb{E}[(Y-\mathbb{E}[Y])^2] = \sum^n_{i=1}(x_i-\mathbb{E}[Y])^2\underbrace{P(Y=x_i)}_{1/n}=\sum^n_{i=1}\frac{(x_i-\bar{x})^2}{n}=\frac{(n-1)}{n}s^2\]

• Bootstrap approximation of sample mean variance:
  • This approaches the true sample mean variance as \(n\rightarrow\infty\)

\[ Var\left(\frac{1}{n}\sum_iY_i\right) = \frac{n-1}{n}\frac{s^2}{n}\rightarrow \frac{s^2}{n}\]

• The following is an R example for the MSE of the sample median:

tit <- read.table('titanic.txt', header = T)
tit <- na.omit(tit)
res <- rep(0,times=1000)
for (i in 1:1000){
  res[i] <- median(sample(tit$Age, replace = T))
}
mean((res-median(tit$Age))^2)

## [1] 0.56975

plot(density(res), type = 'l')