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)
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')
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)
(1-thetacan)^(b+N-x)-1
is missing a parenthesis in the acceptance probability and so does(1-theta)^(b+N-x)-1
$\endgroup$