Problem 1

Generate data from a Cauchy random variable with density \(f(x)=1/[\pi(1+x^2)]\) using two different method discussed in class. Verify your solution by drawing the histogram with the density.

  1. the general inverse method

\[ \begin{aligned} &F_X(x) =\int_{-\infty}^x 1/[\pi(1+t^2)] dt =\\ &\frac1\pi\arctan t|_{-\infty}^x = (\arctan x +\pi/2)/\pi =\frac12+\frac1\pi\arctan x \\ &y=\frac12+\frac1\pi\arctan x \\ &\pi(y-1/2)=\arctan x \\ &x= \tan [\pi(y-1/2)] \end{aligned} \]

The Cauchy rv give occasionally very large values, so I will draw it from -5 to 5:

df=function(x) dcauchy(x)/diff(pcauchy(c(-5, 5)))
B=1e4
x=tan(pi*runif(B,-0.5, 0.5))
x=x[x>-5&x<5]
hist(x, 50, main="", freq=FALSE)
curve(df, -5, 5, add=TRUE, col="blue", lwd=2)

  1. the Metropolis-Hastings algorithm

A a proposal I will use a simple uniform neighbor, so

\[\frac{b_{x}q_{xx_n}}{b_{x_n}q_{x_nx}}=\frac{1+x_n^2}{1+x^2}\]

finalp1MH <- function(B, eps) {
  x=rep(0, B)
  for(i in 2:B) {
      y=x[i-1]+runif(1, -eps, eps)
      if(runif(1)<(1+x[i-1]^2)/(1+y^2)) x[i]=y
      else x[i]=x[i-1]
  }
  x
}

As the next calculaion shows, this takes a while to reach the stationary distribution:

B=1e5
x=matrix(0, B, 5)
for(i in 1:5) x[, i]=finalp1MH(B, eps=5)
x1=apply(x, 2, cumsum)/1:B
plot(1:B, x1[, 1], pch=".", ylim=range(x1))
for(i in 2:5) lines(1:B, x1[, i])

so we will just use the last 10000 of each:

x=c(x[(B-10000):B, ])
x=x[x>-5&x<5]
hist(x, 50, main="", freq=FALSE)
curve(df, -5, 5, add=TRUE, col="blue", lwd=2)

Problem 2

Generate data from the density proportional to

\[g(x, y)=\sin(\pi x/2)\frac{\Gamma(2x+2)}{\Gamma(x+1)^2}[y(1-y)]^x;0<x,y<1\]

using the runif command only. Verify your solution by drawing the histogram of the marginals with their densities.

\[ \begin{aligned} &f_X(x) =c\int_0^1 \sin(\pi x/2)\frac{\Gamma(2x+2)}{\Gamma(x+1)^2}[y(1-y)]^x dy = \\ &c\sin(\pi x/2) \int_0^1 \frac{\Gamma(2(x+1))}{\Gamma(x+1)^2}[y(1-y)]^{(x+1)-1} dy = \\ &c\sin(\pi x/2) \\ \end{aligned} \]

because the integrand is the density of a Beta(x+1,x+1)

\[ \begin{aligned} &f_x(x) = c\sin(\pi x/2); 0<1<x\\ &F_x(x) = \int_0^x c\sin(\pi t/2) dt = \\ &-\frac{2c}{\pi} \cos(\pi t/2)|_0^x = \frac{2c}{\pi}[1-\cos(\pi x/2)]\\ &1=F_x(1)=\frac{2c}{\pi}[1-\cos(\pi/2)]=\frac{2c}{\pi}\\ &y=1-\cos(\pi x/2) \\ &x=2\arccos(1-y)/\pi \end{aligned} \]

so we can generate X using the inverse method.

Now for Y:

\[ \begin{aligned} &f_{Y|X=x} =\frac{f(x,y)}{f_x(x)} = \\ &\frac{\frac\pi{2}\sin(\pi x/2)\frac{\Gamma(2x+2)}{\Gamma(x+1)^2}[y(1-y)]^x}{\frac\pi{2}\sin(\pi x/2)} = \\ &\frac{\Gamma(2x+2)}{\Gamma(x+1)^2}[y(1-y)]^x \end{aligned} \]

For this we can use accept-reject. Note the Beta(a,a) density has its maximum at 0.5, so

rfy.x <- function(x) {
  M=dbeta(0.5, x+1, x+1)
  repeat {
    y=runif(1)
    if(runif(1)<dbeta(y, x+1, x+1)/M) break
  }
  y
}

so now

B=1e4
x=2*acos(runif(B))/pi
y=0*x
for(i in 1:B) y[i]=rfy.x(x[i])

Let’s check:

hist(x, 50, freq=FALSE, main="")
curve(pi/2*sin(pi*x/2), 0, 1, add=TRUE, col="blue", lwd=2)

The marginal density of y we need to find numerically

fxy <- function(x, y) {
  pi/2*sin(pi*x/2)*dbeta(y, x+1, x+1)
}
integrate(fxy, 0, 1, x=1)$value
## [1] 1.570796
hist(y, 50, freq=FALSE, main="")
t=seq(0, 1, length=250)
z=0*t
for(i in 1:250) z[i]=integrate(fxy, 0, 1, y=t[i])$value
lines(t, z,  add=TRUE, col="blue", lwd=2)

Problem 3

Use simulation to find

\[\int_{-\infty}^\infty \int_{x-1}^{x+1} \sqrt{x^2+y^2} e^{-x^2} \text{ } dy \text{ }dx\]

\[ \begin{aligned} &\int_{-\infty}^\infty \int_{x-1}^{x+1} \sqrt{x^2+y^2}e^{-x^2} dy dx = \\ &\int_{-\infty}^\infty \int_{-\infty}^\infty \sqrt{x^2+y^2} I_{[x-1, x+1]}(y) e^{-(\sqrt 2 x)^2/2} dy dx = \\ &2\sqrt{\pi}\int_{-\infty}^\infty \int_{-\infty}^\infty \sqrt{x^2+y^2} \frac12 I_{[x-1, x+1]}(y)\frac1{\sqrt{2\pi\times0.5}} e^{-\frac{x^2 }{2\times 0.5} } dy dx = \\ &2\sqrt{\pi}E[\sqrt{X^2+Y^2}] \end{aligned} \]

where \(X\sim N(0, \sqrt{0.5})\) and \(Y|X=x\sim U[x-1,x+1]\), so

B=1e4
x=rnorm(B, 0, sqrt(1/2))
y=runif(B,x-1, x+1)
2*sqrt(pi)*mean(sqrt(x^2+y^2))
## [1] 3.488426

Numerical check:

f <- function(x, y) sqrt(x^2+y^2)*exp(-x^2)*ifelse(abs(y-x)<1, 1, 0)
double.integral(f, -c(Inf, Inf), c(Inf, Inf))
## [1] 3.495948