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

    # -----------------------------------------
    # Treatment cost example - Example 3.5 from 
    # Brazzale, Davison and Reid (2007)
    # -----------------------------------------
y = c(30, 172, 210, 212, 335, 489, 651, 1263, 1294, 1875,
      2213, 2998, 4935, 
      121, 172, 201, 214, 228, 261, 278, 279, 351, 561,
      622, 694, 848, 853, 1086, 1110, 1243, 2543)
x = c(rep(1,13), rep(2,18))
y = y / 1000

s1 = sum(y[x==1])
s2 = sum(y[x==2])
n=13
m = 18

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


#================================
 loglik <- function( param, y, x)
#================================
{
    # -----------------------------------
    # Log-likelihood for exponential with
    # ratio of means
    # -----------------------------------

    mu1 = param[1]
    mu2 = param[2]
    -( -13*log(mu1) - 18*log(mu2) - 
                 sum(y[x==1])/mu1 - sum(y[x==2])/mu2)
}


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

   param[1]/param[2]
}

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

   param[1] - param[2]
}



param = c(s1/n, s2/m)
truese = matrix( c(0.72050210941595216409, 0.95566632345012524281,
                   0.95566632345012524281, 0.38720003344196362612), 2, 2)

der1 <- g2(transfm1, param)
der2 <- g2(transfm2, param)
dermat = matrix(c(der1,der2), 2, 2, byrow=TRUE)
corr <- dermat %*% solve(h1(loglik,param,y,x)) %*% t(dermat)
tmp = (sqrt(diag(corr)))
err1 = log10( max( abs( tmp - diag(truese)) / diag(truese) ))

der1 <- g2(transfm1, param)
der2 <- g2(transfm2, param)
dermat = matrix(c(der1,der2), 2, 2, byrow=TRUE)
corr <- dermat %*% solve(h2(loglik,param,y,x)) %*% t(dermat)
tmp = (sqrt(diag(corr)))
err2 = log10( max( abs( tmp - diag(truese)) / diag(truese) ))

der1 <- g2(transfm1, param)
der2 <- g2(transfm2, param)
dermat = matrix(c(der1,der2), 2, 2, byrow=TRUE)
corr <- dermat %*% solve(h3(loglik,param,y,x)) %*% t(dermat)
tmp = (sqrt(diag(corr)))
err3 = log10( max( abs( tmp - diag(truese)) / diag(truese) ))

der1 <- g4(transfm1, param)
der2 <- g4(transfm2, param)
dermat = matrix(c(der1,der2), 2, 2, byrow=TRUE)
corr <- dermat %*% solve(h4(loglik,param,y,x)) %*% t(dermat)
tmp = (sqrt(diag(corr)))
err4 = log10( max( abs( tmp - diag(truese)) / diag(truese) ))

der1 <- g2(transfm1, param)
der2 <- g2(transfm2, param)
dermat = matrix(c(der1,der2), 2, 2, byrow=TRUE)
corr <- dermat %*% solve(hessian(loglik,param,,,y,x)) %*% t(dermat)
tmp = (sqrt(diag(corr)))
err5 = log10( max( abs( tmp - diag(truese)) / diag(truese) ))

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