Problem 1

Generate rv’s from a standard normal using the slice sampler.

\[ \exp(-x^2/2)=u \]

\[ \begin{aligned} &\exp(-x^2/2)=u \\ &-x^2/2=\log(u) \\ &x= \pm \sqrt{-2\log(u)} \\ \end{aligned} \]

hw9p1 <- function(B=1e4) {
  f <- function(x) exp(-x^2/2)
  x <- rep(0, B)
  for(i in 2:B) {
    u <- runif(1, 0, f(x[i-1]))
    x[i] <- runif(1, -sqrt(-2*log(u)),
                      sqrt(-2*log(u)))
  }
  x
}
x <- hw9p1()
hist(x, 50, freq=FALSE, main="")
curve(dnorm(x), -3, 3, col="blue", lwd=2, add = TRUE)

Problem 2

Generate rv’s from a density \(f(x,y)=xy+1,0<x,y<1\) using the slice sampler.

Say we have u, then the problem is to sample uniformly form the set of (x,y) with \(xy+1>u\). But what is this set?

\[ \begin{aligned} &xy+1 > u\\ &y > (u-1)/x \end{aligned} \] so we generate such a point by accept-reject. Now

hw9p2 <- function(B=1e4) {
  f <- function(x) prod(x)+1
  x <- matrix(0, B, 2)
  x[1, ] <- c(1/2, 1/2)
  for(i in 2:B) {
    u <- runif(1, 0, f(x[i-1, ]))
    repeat {
      xy <- runif(2)
      if(xy[2]>(u-1)/xy[1]) break
    }
    x[i, ] <- xy
  }
  x
}
x <- hw9p2()

\[ \begin{aligned} &f(x) = \int_0^1 xy+1dy=xy^2/2+y|_0^1 = x/2+1\\ &\int_0^1 x/2+1dx = x^2/4+x|_0^1 = 1/4+1 = 5/4 \end{aligned} \]

par(mfrow=c(1, 2))
hist(x[, 1], 50, freq=FALSE, main="")
curve(4/5*(x/2+1), 0, 1, col="blue", lwd=2, add = TRUE)
hist(x[, 2], 50, freq=FALSE, main="")
curve(4/5*(x/2+1), 0, 1, col="blue", lwd=2, add = TRUE)

Problem 3

Below is a sample from a Geometric random variable with r, that is \(P(X=k)=r(1-r)^{k-1};k=1,2,..\). Let’s say we know that \(0.2<r<0.4\), so we will do a Bayesian analysis using a prior \(\pi(r) = U[0.2, 0.4]\). Find a 90% credible interval for r using the quantiles of data simulated from the posterior distribution.

## x
##  1  2  3  4  5  6  7  8 14 
## 15 11  5  8  5  2  1  2  1

We have

\[ \begin{aligned} &\pi(r) = 2 I_{[0.2,0.4]}(r) \\ &f(x_1,..,x_n|r) =\prod_{i=1}^n \left[ r(1-r)^{x_i-1}\right] = r^n(1-r)^{\sum x_i-n}\\ &f(x_1,..,x_n,r) = r^n(1-r)^{\sum x_i-n} 2 I_{[0.2,0.4]}(r)\\ &f(x_i|x_1,.,x_{i-1}, x_{i+1},..,x_n,r) \propto (1-r)^{x_i-1} \\ &f(r|x_1,..,x_n)\propto Beta(n+1, \sum x_i-n+1) I_{[0.2,0.4]}(r)\\ \end{aligned} \]

so all the marginals are easy to generate, so

hw11p2=function(x, B=1e4) {
  
  n=length(x)
  r=rep(1/mean(x), B)
  for(i in 2:B) {
    x=rgeom(n, r[i])+1
    repeat {
      r[i] = rbeta(1, n+1, sum(x)-n+1)
      if(r[i]>0.2 & r[i]<0.4) break 
    }  
  }
  r
}
x=rep(c(1,  2,  3,  4,  5,  6,  7,  8, 14), 
      c(15, 11,  5,  8,  5,  2,  1,  2,  1))
r=hw11p2(x)
hist(r, 100)

round(quantile(r, c(0.05, 0.95)), 3)
##    5%   95% 
## 0.244 0.387