ECON546 Lab9 Gibbs Sampler

Gibbs sampling, in its basic incarnation, is a special case of the Metropolis-Hastings algorithm. The point of Gibbs sampling is that given a multivariate distribution it is simpler to sample from a conditional distribution than to marginalize by integrating over a joint distribution.

Monte Carlo integration- YouTube

---

MCMC visualization java

Gibbs sampling - YouTube

MCMC and the Metropolis Hastings algorithm- YouTube

Example on slides

MCMC Visualizationi r

---

library(ggplot2)
set.seed (123456)
nrep<- 55000
burn<- 5000
margy<- vector(length=nrep)
margx<- vector(length=nrep)
n<- 10
y<- 0.1      # TRY DIFFERENT INITIAL VALUES FOR Y BY CHANGING THIS NUMBER
a<- 2 # prior shape 1 for beta
b<- 4 # prior shape 2  for beta

(a) What is the conditional distribution for X, given Y?

https://en.wikipedia.org/wiki/Binomial_distribution

\[f(k;n,p)=\Pr(X=k)={\binom {n}{k}}p^{k}(1-p)^{n-k}\]

(b) What is the conditional distribution for Y, given X?

\[f(x; \alpha ,\beta )=\mathrm{constant} \cdot x^{\alpha -1}(1-x)^{\beta -1}\]

$${}$$

https://en.wikipedia.org/wiki/Beta_distribution

---

(c) obtain the marginal distributions of X and Y.

from conditional distribution converge to marginal distribution.

for (i in 1:nrep)  {
x<- rbinom(1, n, y)
margx[i]<- x
y<- rbeta(1, x+a,n-x+b)
margy[i]<- y
}

---

(f) Summarize the features of the marginal distributions for X and Y, and plot the marginal p.m.f. for X and the marginal p.d.f. for Y.

Graph

Posterior distribution

# 
hist(margx[(burn+1):nrep], col="red",prob=TRUE,main="Marginal Mass Function for X", xlab="x", ylab="p(x)")

ggplot()+geom_histogram(mapping = aes(x =margx[(burn+1):nrep] ), fill = "red")+
  labs(y="p(x)", x = "x", title ="Marginal Mass Function for X")

plot(density(margy[(burn+1):nrep]),col="blue",main="Marginal Density for Y", xlab="y", ylab="p(y)")

ggplot()+geom_density(mapping = aes(x = margy[(burn+1):nrep]))+
  labs(y="p(y)", x = "y", title ="Marginal Density Function for y")

---

(d) Establish that E(X) = 3.337, E(Y) = 0.335, Var.(X) = 5.132, and Var.(Y) = 0.032.

summary(margx[(burn+1):nrep])

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   2.000   3.000   3.337   5.000  10.000 

summary(margy[(burn+1):nrep])

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
0.001594 0.193600 0.315300 0.334500 0.457600 0.938400 

var(margx[(burn+1):nrep])

[1] 5.131855

var(margy[(burn+1):nrep])

[1] 0.03214491

---

(e) Show that these results do not depend on the initial value assigned for Y.

y<- 0.1      # TRY DIFFERENT INITIAL VALUES FOR Y BY CHANGING THIS NUMBER

---

(f) Is the Burn-in period long enough?

Rolling mean diagnostics:

rmean_x<- vector(length=burn)
rmean_y<- vector(length=burn)
for (i in 1:burn) {
rmean_x[i]<- mean(margx[1:i])
rmean_y[i]<- mean(margy[1:i])
}

plot(rmean_x, col="red", main="Rolling Means for X",xlab="Replications", ylab="Mean of X")

ggplot()+geom_point(mapping = aes(x = 1:length(rmean_x), y =rmean_x ), col="red")

plot(rmean_y, col="blue", main="Rolling Means for Y", xlab="Burn-in Replications", ylab="Mean of Y")

ggplot()+geom_point(mapping = aes(x = 1:length(rmean_y), y =rmean_y ), col="blue")