## com108c of november 2021,
#  estimating an ensemble of Poisson parameters,
#  via hierarchical setup, for stk 4021 

#  (theta, y) from 
#  gamma(theta, a,b) f(y, theta) = gamma(theta,a+y,b+1) f_marg(y,a,b)
#  with (b^a\over \Gamma(a)) (\Gamma(a+y)\over (b+1)^(a+y)) (1/y!)   

#  need mcmc for (a,b) from the posterior,
#  after which we simply do 
#  theta_i given (data,a,b) is gamma(a+y_i,b+1) 

#  \prior(a,b) \prod_{i=1}^n f_\marg(y_i, a, b) 

#  here i do the mcmc with (a, b) = (c theta0, c),
#  for better chain and better interpretation (and better prior, if needed)

nn <- 40
# yy <- c(rpois(round(0.5*nn),4.44), rpois(nn-round(0.5*nn),8.88))
# yy <- sort(yy) 

yy <- c(0, 2, 2, 3, 3, 3, 3, 3, 3, 3,  
        4, 4, 5, 5, 5, 5, 5, 6, 6, 6,  
        6, 6, 6, 7,  7, 7,  7,  7,  8,  
        8,  9,  9, 10, 10, 11, 12, 12, 12, 15, 18)
 
check <- var(yy)/mean(yy) 
print(c(1+sqrt(2/nn),check))

# the ifelse thing can be used, 
# to avoid logL(a,b) calculations when b < 0: 
# x <- 1.234
# ifelse(x > 0, sqrt(x), -33)

# i take care to construct logL(c(theta0,c)) such that 
# it returns something very long down if a or b is negative 

logL <- function(para)
{
theta0 <- para[1]
cc <- para[2]
# 
a <- cc*theta0
b <- cc 
# 
aux <- ifelse(1*(a > 0)*(b > 0) == 1, 
sum( lgamma(a+yy)-lgamma(a)-(a+yy)*log(b+1)+a*log(b) - lgamma(yy+1) ), -888)
aux 
}

logL(c(10,2))

minuslogL <- function(ab)
{-logL(ab)}

nils <- nlm(minuslogL,c(10,2),hessian=T)
ML <- nils$estimate
Jhat <- nils$hessian
se <- sqrt(diag(solve(Jhat)))
showme <- cbind(ML,se)
print(round(showme,4))

#         ML     se
#[1,] 6.5750 0.5814 theta0 
#[2,] 0.9471 0.4289 c 

sim <- 10^5 # 10^5 ought to be more than enough 
mcmc <- 0*(1:sim) %*% t(1:2) 

# for proposal: 
vv1 <- rnorm(sim)
vv2 <- rnorm(sim) 

uu <- runif(sim)
ok <- 0*(1:sim)
pracc <- 0*(1:sim) 

mcmc[1, ] <- c(ML[1], ML[2]) # starting somewhere 

suitable <- 0.25

for (ss in 2:sim)
{
old <- mcmc[ss-1, ]
prop <- old + c(suitable*se[1]*vv1[ss],suitable*se[2]*vv2[ss]) 
# 
ratio <- exp(logL(prop) - logL(old))
pracc[ss] <- min(c(1,ratio)) * 1*(prop[2] > 0)
ok[ss] <- 1*(uu[ss] <= pracc[ss])
# 
mcmc[ss, ] <- prop*ok[ss] + old*(1-ok[ss])
}

plot(mcmc[ ,1]) 

theta0sim <- mcmc[ ,1]
csim <- mcmc[ ,2]

plot(theta0sim,csim) 

# note: (a, b) = (c theta0, c) 

thetasim <- 0*(1:sim) %*% t(1:nn) 

for (i in 1:nn)
{
thetasim[ ,i] <- rgamma(ss,csim*theta0sim + yy[i],csim+1)
}

quant05 <- 0*(1:nn)
quant50 <- 0*(1:nn)
quant95 <- 0*(1:nn)
means <- 0*(1:nn)
stdevs <- 0*(1:nn) 

for (i in 1:nn) 
{
quant05[i] <- quantile(thetasim[ ,i],0.05)
quant50[i] <- quantile(thetasim[ ,i],0.50)
quant95[i] <- quantile(thetasim[ ,i],0.95)
means[i] <- mean(thetasim[ ,i])
stdevs[i] <- sd(thetasim[ ,i])
}

## making a plot, with raw data (same as ML estimates)
#  plus Bayes estimators and 0.05, 0.50, 0.95 intervals: 

par(bty="l")
par(cex.lab=1.4142) 

plot(1:nn,yy,
   xlab="the studies, ordered in outcome",
   ylab="raw data and intervals")

## matpoints(1:nn,means,pch=2,col=4)
matlines(1:nn,cbind(quant05,quant50,quant95),type="l",
   lwd=3.0,col=c(2,4,2),lty=c(2,1,2))

okthen <- cbind(1:nn,yy,quant05,quant50,quant95)
print(round(okthen,4))