0
$\begingroup$

Let's say that we have calculated the posterior distribution of a parameter of interest given the data of a binomial experiment $N=70,x=34$ which the probability of event occurrence $\theta$ follows the Beta distribution with parameters $\alpha=4.4,\beta =6.6$.So we have a posterior distribution :

$$p(\theta|x) \propto f(x|\theta)f(\theta) \propto \theta^{(\alpha+x)-1}(1-\theta)^{(\beta+N-x)-1}$$

Therefore we know that the un-normalized posterior distribution of parameter $\theta$ given the data follows the Beta distribution, i.e $$\theta\sim Beta (38.4,42.6) $$

For classical statistics (frequentist) one can calculate the probability of occurrence as: $$\theta_{classical} = \frac{X}{N} \approx 0.48$$ and a Bayesian would do : $$\theta_{Bayesian} = \mathbb{E}[\theta|x] = \frac{\alpha+x}{\alpha+\beta+N} = \frac{38.4}{38.4+42.6}=0.4741$$

I am now plotting the prior and the posterior distribution :

alpha = 4.4
beta  =  6.6
n = 70
x = 34
theta <- seq(0,1, by = 0.01) # set up grid for plotting
plot(theta, dbeta(x=theta, alpha, beta), type = 'l', lwd = 2, col = 'orange',
     ylim = c(0, 10), xlab = expression(theta), 
     ylab = expression(paste('p(', theta, '|y)')))
lines(theta, dbeta(theta, alpha + x, beta + n-x), 
      type = 'l', lwd = 2, col = 'violet')
legend('topright', inset = .02, legend = c('prior', 'posterior'),
       col = c('orange', 'violet'), lwd = 2)

enter image description here

So far so good and it seems extremely easy regarding the calculations. Conceptually with I have done is that I used a prior knowledge of this specific experiment in order to calculate the probability of occurrence.I used a conjugate prior for that and the posterior distribution has an analytical form. My first question is :If I use a conjugate prior is it necessary to calculate the integral in the denominator of Bayes Law ?

Bays Law : $$p(\theta|x) =\frac{f(x|\theta)f(\theta)}{\int f(x|\theta)f(\theta)d\theta }$$

The second question is: Can I use the Metropolis Hastings algorithm in order to draw samples from the posterior?

I took the liberty of doing so in R (code below) and I used a starting value of $\theta^{0}$ from the prior distribution $Beta(\alpha=4.4,\beta = 6.6)$ and the proposal distribution of the algorithm to be uniform on $[0,1]$. But the resulting $\theta$ from the convergence plot is way off the 0.47


mh = function(N,x,a,b,nburn=0,ndraw=1000){ 

  #initial value: drawn from prior
  theta = rbeta(1,4.4,6.6)

  # vector of recorded draws:
  draws = numeric(ndraw)

  # counter for acceptance probability:
  accept = 0  

  # MCMC LOOP FOLLOWS:
  it = -nburn
  while(it < ndraw){
    it = it+1;  

  # draw thetacan from Uniform(0,1):
  thetacan = runif(1,0,1)

  # compute acceptance probability:
  alpha = ((thetacan^((a+x)-1))*( (1-thetacan)^(b+N-x)-1 ))/
          ((   theta^((a+x)-1))*( (1-   theta)^(b+N-x)-1 ))

  # draw u ~ Uniform (0,1):
  u = runif(1,0,1)

  # if u<p, take thetacan as next value of chain:
  if (u < alpha){
    accept = accept + 1 
    theta = thetacan}
  
  # after burn-in record theta:
  if(it>0){draws[it] <- theta}
  }
  # END MCMC
  return(draws)
}

ndraw = 10000
theta=mh(N=70,x=34,a=4.4,b=6.6,nburn=0,ndraw=ndraw) 
plot(cumsum(theta)/1:ndraw,type='l')


enter image description here

My last question is what I made wrong in the mcmc code and it is way off the 0.47? (if it is ok to calculate the $\theta$ parameter in this way)

$\endgroup$
5
  • 1
    $\begingroup$ See stats.stackexchange.com/q/129666/7224 $\endgroup$
    – Xi'an
    Commented Dec 11, 2022 at 13:23
  • 1
    $\begingroup$ R error : (1-thetacan)^(b+N-x)-1 is missing a parenthesis in the acceptance probability and so does (1-theta)^(b+N-x)-1 $\endgroup$
    – Xi'an
    Commented Dec 11, 2022 at 13:27
  • 1
    $\begingroup$ "the un-normalized posterior distribution of parameter θ given the data follows the Beta distribution" : no, the posterior distribution IS the Beta distribution and its normalising constant is provided from the fully known density of the Beta. Hence in conjugate settings one rarely need compute the normalising constant. $\endgroup$
    – Xi'an
    Commented Dec 11, 2022 at 13:29
  • 1
    $\begingroup$ @Xi'an thank you for comment professor.So if I do actually implement the MH algorithm is useless in this conjugate setting? $\endgroup$ Commented Dec 11, 2022 at 13:51
  • 1
    $\begingroup$ Since the conjugate posterior is from a standard family of distributions, there should exist an already implemented code to simulate from this family of distributions. In almost all situations, this should prove more efficient than running a generic MH algorithm. $\endgroup$
    – Xi'an
    Commented Dec 11, 2022 at 14:42

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.