#==============
# Preliminaries
#==============

require(MASS)
require(nlme)
require(numDeriv)
source("hessian.txt")


#===================================
 loglik1 <- function( param, y, x) {
#===================================

        # Log-likelihood function for smoking data

#    gamma = 1 / (1 + exp(-param[1]))
    gamma = param[1]
     -(y * log(gamma) + (x-y) * log(1-gamma))
}

#=============================
 transfm <- function( param) {
#=============================

   param[1]/(1-param[1])
}

#====================
# Davison Example 4.6
#====================

    #------------------
    # Mobile phone data
    #------------------
y = 157
n = 181
x = n

truese = 1.4337461049706529258
param = 157/181


hess1 <- h1(loglik1, param, y, x)
corr1 <- hess2corr(hess1)
corr1 <- corr1 * g2(transfm, param)

hess2 <- h2(loglik1, param, y, x)
corr2 <- hess2corr(hess2)
corr2 <- corr2 * g2(transfm, param)

hess3 <- h3(loglik1, param, y, x)
corr3 <- hess2corr(hess3)
corr3 <- corr3 * g2(transfm, param)

hess4 <- h4(loglik1, param, y, x)
corr4 <- hess2corr(hess4)
corr4 <- corr4 * g4(transfm, param)

hess5 <- hessian(loglik1,param,,,y,x)
#hess5 <- hessian(loglik1,param,,method.args=list(eps=1e-4, d=0.1, r=6, v=2),y,x)
corr5 <- hess2corr(hess5)
#corr5 <- corr5 * g2(transfm, param)
corr5 <- corr5 * grad(transfm,param)

err1 = log10(abs( (corr1 - truese) / truese))
err2 = log10(abs( (corr2 - truese) / truese))
err3 = log10(abs( (corr3 - truese) / truese))
err4 = log10(abs( (corr4 - truese) / truese))
err5 = log10(abs( (corr5 - truese) / truese))

cat(err1,err2,err3,err4,err5,"\n")

#tim = c(t1[3],t2[3],t3[3],t4[3],t5[3])
#print(tim)
#print(tim/tim[3])


#stop()