Homework Ten - Simio and Simulation: Chapter 8 - Discussion

Complete problem 5 from Chapter 8 in Kelton.

#create reproducability
set.seed(123)

Discussion

Attempt to implement a sampler for 9.3 or 9.4 in Rizzo’s text.

9.3

Use the Metropolis-Hastings sampler to generate random variables from a standard Cauchy distribution. Discard the first 1000 of the chain, and compare the deciles of the generated observations with the deciles of the standard Cauchy distribution (see qcuchy or qt with df=1). Recall that a Cauchy(\(\theta\), \(\eta\)) distribution has density function

\(f(x) = \frac{1}{\theta\pi(1+[(x - \eta)/\theta]^2)}\), \(-\infty\) < x < \(\infty\), \(\theta\) > 0

The standard Cauchy has the Cauchy(\(\theta\) = 1, \(\eta\) = 0) density. (Note that the standard Cauchy density is equal to the Student t density with one degree of freedom.)

Metropolis-Hastings Sampler generates the Markov chain {X(0), X(1)…} as follows: 1. Choose a proposal distribution g(.|X(t)) (subject to regularity conditions). 2. Generate X(0) from a distribution g. 3. Repeat (until the chain has converged to a stationary distribution according to some criterion): a. Generate Y from g(.|X(t)) b. Generate U from Uniform(0,1) c. if U <= d. Increment t.

#sample size
n <- 10000
x <- numeric(n)
u <- runif(n)
rejects <- 0

for (i in 2:n){
  
  j <- x[i-1]
  y <- rnorm(1)
  
  if (u[i] <= dcauchy(y) * dnorm(j, mean = y)/dcauchy(j) * dnorm(y, mean = j)) {
    x[i] <- y
  }
  else {
    x[i] <- j
  }
}

#remove first thousand results
x1 <- x[1001:length(n)]
observations <- quantile(x, seq(0, 1, 0.1))
expectations <- qcauchy(seq(0, 1, 0.1))
decile <- data.frame(observations, expectations)
decile
##      observations expectations
## 0%   -2.020557058         -Inf
## 10%  -0.647321270   -3.0776835
## 20%  -0.422131687   -1.3763819
## 30%  -0.268546251   -0.7265425
## 40%  -0.118943640   -0.3249197
## 50%  -0.004460459    0.0000000
## 60%   0.111538963    0.3249197
## 70%   0.239866735    0.7265425
## 80%   0.419217723    1.3763819
## 90%   0.703869507    3.0776835
## 100%  1.661924638          Inf

9.4

Implement a random walk Metropolis sampler for generating the standard Laplace distribution. For the increment, simulate from a normal distribution. Compare the chains generated when different variances are used for the proposal distribution. Also, compute the acceptance rates of each chain.