# Posterior simulation under the Normal-Normal hierarchical model
#
# Computation of the posterior distribution of (theta,mu,tau) via 
# simulation following the factorization:
# 
# p(theta,mu,tau|y) propto p(tau|y)p(mu|tau,y)p(theta|mu,tau,y)
#
# 1) Computation of density p(tau|y) as in (5.21) with 
# p(tau) propto 1, at p. 139, Gelman's book 
#
# Arguments: e=c(y1,...,y_n); v=c(sigma2_1,...,sigma2_n);
# x.tau=seq(0,1,.001).
#
# NOTE  In beta-blockers example y_j values are the empirical logits 
# and sigma2_j values are the approximate sampling variances, 
# computed as from (5.23) and (5.24), p. 150, Gelman's book.
#
dtau.r_function(e,v,x.tau)
{
n_length(e)
d.tau_c()
for(i in 1:length(x.tau))
{
prec_1/(x.tau[i]^2+v)
v.mu_( sum( prec ) )^(-1)
mu.bar_(prec%*%e) * v.mu

d.tau[i]_sqrt(v.mu)

for(j in 1:n)
{ d.tau[i]_d.tau[i]*
  (v[j]+x.tau[i]^2)^(-1/2)*
  exp(-((e[j]-mu.bar)^2)/(2*(v[j]+x.tau[i]^2))) 
}}
d.tau
}
#
# Draw tau from p(tau|y)
# Inverse cdf method 
# Arguments: d.tau=dtau.r output; m=#draws
#
tau.r_function(d.tau,m)
{
cn_sum(d.tau)
dn.tau_d.tau/cn
f.tau_cumsum(dn.tau)

tau_c()
for(i in 1:m)
{
u_max(runif(1),f.tau[1])
cat(u,"\n")
tau[i]_x.tau[length(f.tau[f.tau<u | f.tau==u])]
}
tau
}
#
# Draw mu from p(mu|tau,y), with (mu|tau) propto 1
# mu|tau,y ~ N(mu.bar,v.mu), as in (5.20), p.139, Gelman's book 
# Arguments: e,v as above; tau=tau.r output
#
mu.r_function(e,v,tau)
{
mu_c()
for(i in 1:length(tau))
{
prec_1/(tau[i]^2+v)
v.mu_( sum( prec ) )^(-1)
mu.bar_(prec%*%e) * v.mu
mu[i]_rnorm(1,mu.bar,sqrt(v.mu))
}
mu
}
#
# Draw from p(theta|mu,tau,y)
# theta|mu,tau,y ~ N(theta_j.bar,v_j.theta), asin (5.17), p.138, 
# Gelman's book 
# Arguments: e,v,mu,tau,m as above
#
theta.r_function(e,v,mu,tau,m)
{
n_length(e)
theta_matrix(,m,n)
for(i in 1:m){
 for(j in 1:n){
  theta.bar_(1/v[j] * e[j] + 1/tau[i]^2 * mu[i])/
            (1/v[j] + 1/tau[i]^2)
  v.theta_(1/v[j] + 1/tau[i]^2)^(-1)
  theta[i,j]_rnorm(1,theta.bar,v.theta)
}}
theta
}