Problem 1

a

Here \(\bar{\theta_1} = 8\) and \(s_{\theta_1} = 4\). The standard deviation due to simulation variability is \(s_{\theta_1} / \sqrt{S} = 4 / \sqrt{100} = .4\)

b

To reduce the standard deviation by \(1/4\), the square root of sample size should be increased by factor \(4\), and the sample size should be increased by factor \(16\), so \(1600\) draws.

c

The asymptotic variance of quantile estimator is \(\frac{p(1-p)}{nf(\hat{q})^2}\). Since it’s normal distribution, the standard error of estimate for \(.025\) and \(.975\) quantiles are equal. Both of them can be calculated as:

q.hat <- qnorm(.025, mean = 8, sd = 4)
sqrt(.025*.975/100/(dnorm(q.hat, 8, 4))^2)
## [1] 1.068524

Problem 2

a

Use uniform prior on \(p_1\) and \(p_2\), the posterior are both Beta densities with parameter \((7, 5)\) and \((11, 11)\).

p1 <- rbeta(1000, 7, 5) 
p2 <- rbeta(1000, 11, 11)
quantile(p1 - p2, probs = c(.025, .975))
##       2.5%      97.5% 
## -0.2853458  0.4173117

There are 680 out of \(1000\) observations that \(p1>p2\), so the probability of \(p1 > p2 = 0.68\).

b

The probability of \(p_1 > p_2\) is

\[ Pr[p_1 > p_2] = \int_0^1 p(p_1\mid y_1) \int_0^{p_1} p(p_2\mid y_2) d p_2 d p_1 \] Where \(p(p_i\mid y_i), i = 1, 2\) are beta density functions.

joint.density <- function(p1) {
  dbeta(p1, 7, 5) * integrate(dbeta, 0, p1, shape1 = 11, shape2 = 11)$value
}
integrate(Vectorize(joint.density), 0, 1)
## 0.6863834 with absolute error < 2.4e-11

Problem 3

a

Since in this Beta mixture, x is restricted within \((0, 1)\), a uniform \(g(x)\) is chosen. The maximum density of the mixture is approximately 1.595.

fx <- function(x) .7 * dbeta(x, 4, 2) + .3 * dbeta(x, 1.5, 3)
gx <- function(x) dunif(x, 0, 1)
optimise(fx, c(0, 1), maximum = TRUE)
## $maximum
## [1] 0.7210772
## 
## $objective
## [1] 1.594117
theta <- runif(1000)
u <- runif(1000, 0, 1.595)
x <- theta[u < fx(theta)/gx(theta)]
length(x)
## [1] 627

b

A Beta(\(c\theta, c(1-\theta)\)) density is chosen as proposal density. Its mean is \(\theta\) and variance is \(\theta(1-\theta)/(c+1)\). Parameter \(c\) is used to control the variance.

dprop <- function(x, theta, c = 3) dbeta(x, c*theta, c*(1-theta))
rprop <- function(theta, c = 3) rbeta(1, c*theta, c*(1-theta))

theta <- vector()
acpt <- 0
theta[1] <- .5
for (i in 2:1500) {
  theta.prop <- rprop(theta[i-1])
  log.r <- log(fx(theta.prop)) - log(fx(theta[i-1])) + 
    log(dprop(theta[i-1], theta.prop)) - log(dprop(theta.prop, theta[i-1]))
  
  if (runif(1) < exp(log.r)) {
    theta[i] <- theta.prop
    if (i > 500) acpt <- acpt + 1
  } else {
    theta[i] <- theta[i-1]
  }
}
theta <- theta[-(1:500)]
acpt/1000
## [1] 0.517

Problem 4

The two conditional densities are both Normal. Their mean and variance are: \[ x_1\mid x_2 \sim N(1 + 2 \cdot 1/2 \cdot (x_2 - 1), (1-1/4) \cdot 4 ) = N(x_2, 3)\] \[ x_2\mid x_1 \sim N(1 + 1/2 \cdot 1/2 \cdot (x_1 - 1), (1-1/4) \cdot 1) = N((x_1 + 3)/4, 3/4) \]

x <- data.frame(x1 = 1, x2 = 1)
for (i in 2:1200) {
  x[i, 1] <- rnorm(1, x[i-1, 2], sqrt(3))
  x[i, 2] <- rnorm(1, (x[i, 1] + 3)/4, sqrt(3/4))
}
x <- x[-(1:200), ]
cov(x)
##           x1        x2
## x1 3.6815272 0.9060638
## x2 0.9060638 0.9650390