############################################
# STAT 740 R examples
# Tim Hanson
# U. South Carolina
# September 22, 2017
############################################

############################################
# compute E(x^4) for x ~ N(1,4)
############################################

f=Vectorize(function(x){x^4*exp(-0.5*(x-1)^2/4)/sqrt(2*pi*4)})
integrate(f,-13,15)

reimann=function(f,a,b,J){
 delta=(b-a)/J
 delta*sum(f(a+delta*(1:J-0.5)))
}

# note that "effective range" more important than J!
reimann(f,-5,7,10)
reimann(f,-5,7,100)
reimann(f,-13,15,10)
reimann(f,-13,15,100)

library(fastGHQuad)
rule=gaussHermiteData(10)
rule

sum((1+rule$x*sqrt(2)*2)^4*rule$w)/sqrt(pi)

############################################
# Ache hunting data
############################################

library(numDeriv)
library(fastGHQuad)
library(lme4) # has glmer
rule=gaussHermiteData(100)
# rule=gaussHermiteData(1) is simple Laplace approximation

d=read.table("http://people.stat.sc.edu/hansont/stat740/ache.txt",header=T)
attach(d) # makes 'age', 'days', and 'monkeys' available to R
n=length(age)
X=cbind(rep(1,n),age)
id=1:47 

m=monkeys; a=age; t=days
# Here's the easy way...
f=glmer(monkeys~age+(1|id),family="poisson",offset=log(days),nAGQ=100)
summary(f)

# Now by hand...
f.glm=glm(monkeys~age,family="poisson",offset=log(days))
theta.old=c(f.glm$coef,0) # fit of GLM provides beta estimate
u=rep(0,47) # obtain Empirical Bayes (EB) estimates of u_i
sigma=1 # need reasonable starting value for sigma
for(i in 1:n){
 lp=X[i,]%*%theta.old[1:2] # use GLM estimate to start
 uo=0
 repeat{ # Newton-Raphson for each u_i
  un=uo-(-t[i]*exp(lp+uo)+m[i]-uo/sigma^2)/(-t[i]*exp(lp+uo)-1/sigma^2)
  if(abs(un-uo)<0.000001){break}
  uo=un
 }
 u[i]=uo
}
theta.old[3]=log(sd(u)) # estimate of log sigma from EB u_i

ll=function(theta){
 y=rep(0,n); sigma=exp(theta[3])
 for(i in 1:n){
  lp=(X[i,]%*%theta[1:2])[1,1] # ith likelihood portion conditional on u_i
  g=function(u){dpois(m[i],t[i]*exp(lp+u))*dnorm(u,0,exp(theta[3]))}
  uo=0
  repeat{ # find mode of ith likelihood portion as function of u
   un=uo-(-t[i]*exp(lp+uo)+m[i]-uo/sigma^2)/(-t[i]*exp(lp+uo)-1/sigma^2)
   if(abs(un-uo)<0.000001){break}
   uo=un
  }
  scale=sqrt(-1/(-t[i]*exp(lp+uo)-1/sigma^2))
  mode=uo
  y[i]=aghQuad(g,mode,scale,rule)
 }
 sum(log(y))
}

it=0
repeat{
 theta.new=theta.old-solve(hessian(ll,theta.old))%*%t(jacobian(ll,theta.old))
 it=it+1; cat("iteration=",it," theta=",theta.new,"\n")
 if(max(abs(theta.new-theta.old))<0.000001){break}
 theta.old=theta.new
}

# posterior standard deviations (not standard errors!)
sqrt(diag(solve(-hessian(ll,theta.new))))

# use of optim
theta.old=c(f.glm$coef,1)
f=optim(theta.old,ll,control=list(fnscale=-1),hessian=T)
cov=sqrt(diag(solve(-f$hessian)))
est=f$par

#######################################
# verifying that the integration works!
#######################################

theta=c(-5.56013,0.04594,log(1.536)) # from glmer

i=3

g=Vectorize(function(u){dpois(m[i],t[i]*exp(X[i,]%*%theta[1:2]+u))
                        *dnorm(u,0,exp(theta[3]))})
x=seq(-7,7,length=1000)
plot(x,g(x),type="l")
integrate(function(x){g(x)},-1,4)

uo=0 # Newton-Raphson to find the mode
sigma=exp(theta[3])
lp=(X[i,]%*%theta[1:2])[1,1]
repeat{
 un=uo-(-t[i]*exp(lp+uo)+m[i]-uo/sigma^2)/(-t[i]*exp(lp+uo)-1/sigma^2)
 if(abs(un-uo)<0.000001){break}
 uo=un
}
scale=sqrt(-1/(-t[i]*exp(lp+uo)-1/sigma^2))
mode=uo
aghQuad(g(x),mode,scale,rule)