Lab Homework Assignment

Code

n=15
a=3
b=5
N=10^5
x=numeric(N)
x[1]=sample(c(0:n),1)
for(i in 2:N){
  y = rbeta(1,x[i-1]+a , n- x[i-1]+b)
  x[i] = rbinom(1,n,y)
  rm(y)
}

px <- function(x, n, a, b){
  choose(n, x)*gamma(a+b)*gamma(x+a)*gamma(n - x + b)/(gamma(a)*gamma(b)*gamma(a + b + n))
}

z <- c(0 : n)

Plot

hist(x, breaks = seq(-0.5, n + 0.5,by = 1),freq = FALSE, col = "yellow", main = "Beta-Binomial distribution",ylab = "P(x)")
points(z, px(z, n, a, b),type = "h", col = 4, lw = 5)
legend("topright", c("MCMC","Theorical p.m.f"), col = c(1,4), lty = c(1, 1))