Consider the random vector (X,Y,Z) with joint density proportional to

\[g(x,y,z)=4xz+6xy^2+8yz^3\]

when \(0<x,y,z<1\), 0 otherwise. Generate data from this rv

  1. using accept-reject
  2. using Metropolis-Hastings
  3. using the Gibbs sampler

Verify that your methods work by drawing the marginals.

\[ \begin{aligned} &f_{XY}(x,y) =\int_0^1 4xz+6xy^2+8yz^3dz=\\ &2xz^2+6xy^2z+2yz^4 |_0^1 = 2x+6xy^2+2y\\ &f_X(x) =\int_0^1 2x+6xy^2+2y dy = \\ &2xy+2xy^3+y^2|_0^1 = 4x+1\\ &\int_0^1 4x+1 dx = 2x^2+x|_0^1=3\\ &f_Y(y) =\int_0^1 2x+6xy^2+2y dx = \\ &x^2+3x^2y^2+2xy|_0^1=3y^2+2y+1\\ &f_{YZ}(y,z) =\int_0^1 4xz+6xy^2+8yz^3dx=\\ &2x^2z+3x^2y^2+8yz^3x|_0^1 = 2z+3y^2+8yz^3\\ &f_Z(z)=\int_0^1 2z+3y^2+8yz^3dy=\\ &2zy+y^3+4y^2z^3|_0^1=4z^3+2z+1 \end{aligned} \]

  1. accept reject

use \((U,V,W)\sim U[0,1]\), then

\[ \max\{\frac{f(x,y,z)}{g(x,y,z)}\} = \max\{(4xz+6xy^2+8yz^3)/3\} =6 \]

hw10.ac <- function(B=1e4) {
  f.cg <- function(x) {
    (4*x[1]*x[3]+6*x[1]*x[2]^2+8*x[2]*x[3]^3)/18
  }
  x <- matrix(0, B, 3)
  for(i in 1:B) {
    repeat {
      u <- runif(3)
      if(runif(1)<f.cg(u)) {
        x[i, ] <- u
        break
      }
    }
  }
  x
}
show.marginals <- function(x) {
   par(mfrow=c(2, 2))  
  hist(x[, 1], 50, freq=FALSE, main="", xlab="x", ylab="")
  curve((4*x+1)/3, 0, lwd=2, col="blue", add=TRUE)
  hist(x[, 2], 50, freq=FALSE, main="", xlab="y", ylab="")
  curve((3*x^2+2*x+1)/3, 0, lwd=2, col="blue", add=TRUE)
    hist(x[, 3], 50, freq=FALSE, main="", xlab="z", ylab="")
  curve((4*x^3+2*x+1)/3, 0, lwd=2, col="blue", add=TRUE)
}
show.marginals(hw10.ac())

  1. Metropolis-Hastings

use the same (U,V,W)!

hw10.MH <- function(B=1e4) {
  f <- function(x) {
    (4*x[1]*x[3]+6*x[1]*x[2]^2+8*x[2]*x[3]^3)/3
  }
  x <- matrix(0, 2*B, 3)
  x[1, ] <- rep(1/2, 3)
  for(i in 2:(2*B)) {
      u <- runif(3)
      if(runif(1)<f(u)/f(x[i-1, ])) 
        x[i, ] <- u
      else
        x[i, ] <- x[i-1, ]
  }
  plot(1:(2*B), cumsum(x[, 1])/1:(2*B), 
       ylim=c(0,1), type="l", ylab="", xlab="Index")
  lines(1:(2*B), cumsum(x[, 2])/1:(2*B))
  lines(1:(2*B), cumsum(x[, 3])/1:(2*B))
  x[-c(1:B), ]
}
show.marginals(hw10.MH())

  1. Gibbs sampler

here we need the conditional distributions:

\[ \begin{aligned} &f_{X|Y=y,Z=z} (x|y,z) =\frac{4xz+6xy^2+8yz^3}{2z+3y^2+8yz^3} \\ &f_{Z|X=x,Y=} (z|x,y) =\frac{4xz+6xy^2+8yz^3}{2x+6xy^2+2y} \\ &f_{XZ}(y,z)=\int_0^1 4xz+6xy^2+8yz^3dy=\\ &4xzy+2xy^3+4y^2z^3|_0^1 = 4xz+2x+4z^3 \\ &f_{Y|X=x,Z=z} (y|x,z) =\frac{4xz+6xy^2+8yz^3}{4xz+2x+4z^3} \end{aligned} \]

Now each of these in turn can be generated using accept-reject:

gen.marg <- function(f, x, y, z) {
  x <- seq(0, 1, length=1000)
  M <- max(f(x, y, z))
  repeat {
    u <- runif(1)
    if(runif(1)<f(u, y, z)/M) return(u)
  }
  
}
fxyz <- function(x, y, z) 4*x*z+6*x*y^2+8*y*z^3
f <- list(fx = fxyz, 
          fy = function(y, x, z) fxyz(x, y, z), 
          fz = function(z, x, y) fxyz(x, y, z))
par(mfrow=c(2, 2))
for(k in 1:3) {
  x <- rep(0, 1000)
  for(i in 1:1000) 
    x[i] <- gen.marg(f[[k]], x=0.5, y=0.5, z=0.5)
  hist(x, 50, freq=FALSE, main="")
}

and now the Gibbs sampler:

hw10.Gibbs <- function(B=1e4) {
  x <- matrix(0, 2*B, 3)
  x[1, ] <- rep(1/2, 3)
  for(i in 2:(2*B)) {
    x[i, 1] <- gen.marg(f[[1]], 
                x=x[i-1, 1], y=x[i-1, 2], z=x[i-1, 3])
    x[i, 2] <- gen.marg(f[[2]], 
                x=x[i, 1], y=x[i-1, 2], z=x[i-1, 3])
    x[i, 3] <- gen.marg(f[[3]], 
                x=x[i, 1], y=x[i, 2], z=x[i-1, 3])    
  }
  plot(1:(2*B), cumsum(x[, 1])/1:(2*B), 
       ylim=c(0,1), type="l", ylab="", xlab="Index")
  lines(1:(2*B), cumsum(x[, 2])/1:(2*B))
  lines(1:(2*B), cumsum(x[, 3])/1:(2*B))
  x[-c(1:B), ]
}
show.marginals(hw10.Gibbs())