This is from task 5.3. in Bayesian Data Analysis 3rd edition by Gelman et al. It deals with a hierarchical model where I am supposed to simulate the posterior distribution for $\tau$, which is the standard deviation for $\theta$. The details on the model is not really important here. My question deals with the normalization of the simulated posterior distribution:
# 5.3
# data
y = c(28,8,-3,7,-1,1,18,12)
s = c(15,10,16,11,9,11,10,18)
J = length(y)
get the likelihood for each value of tau
p.tau = function(tau,y,s) {
num = sum(1/(s^2+tau^2))
mu.hat = sum(1/(s^2+tau^2)y / num)
V.mu = 1 / num
p = sqrt(V.mu)prod(1/sqrt(s^2+tau^2))exp(-sum((y-mu.hat)^2 / (2(s^2+tau^2))))
}
calculate marginal distribution for tau
tau = seq(0.001,30,length=100)
p = rep(NA, length(tau))
for (i in 1:length(tau)) {
p[i] = p.tau(tau[i],y,s)
}
delta = diff(tau[1:2])
p = p/sum(p*delta)
plot(tau,p,type="l")
The marginal distribution needs to be normalized, and I know a density is supposed to sum to 1. This would normally be done by
p = p/sum(p)
However, here we were asked to do
p = p/sum(p*delta)
with
delta = diff(tau[1:2])
This does not make the density sum to 1. So why are we supposed to do it this way?
