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

    # ------------------------------------------
    # Radioimmunoassay example - Example 5.4 from 
    # Brazzale, Davison and Reid (2007). 
    # ------------------------------------------
ria.data <- read.table("ria.txt", header=TRUE)
x = ria.data$conc
y = ria.data$count

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


#===============================
 loglik <- function( beta, y, x)
#===============================
{
    # -------------------------------------------------
    # Log-likelihood function for radioimmunoassay data
    # -------------------------------------------------

    mu = beta[1] + (beta[2]-beta[1]) / (1 + (x/beta[4])^beta[3])
    wt = mu^beta[5]
    loglik = sum( (y-mu)*(y-mu) / (2 * (beta[6]) * wt) + 
                0.5 * (log(wt) + log(beta[6])))
}

param = c(1.802, 24.575, 1.899, 334.956, 2.095, exp(-8.094))

    # ---------------------------------
    # True obs information matrix from 
    # high-precision MAPLE calculations
    # ---------------------------------
truese = matrix( c(
.196816874386331274239769429599, -.402513119563683809366030206156, 
.909976600860891394246996163811, -.192123665432045178897238849952, 
.180009326296545764367708277240, -.175103699996695180664621673847, 
-.402513119563683809366030206156, .205493574761820455519312642706,
 -.535513263932020245233646582011, -.540265562572632296551266978691, 
-0.948059313982808533705260942973e-1, 0.874409142585676611783165499024e-1,
.909976600860891394246996163811, -.535513263932020245233646582011, 
0.741162389758848985064111199265e-1, .143620150765492549865289851418, 
.175875688555801634799614555557, -.170820256565046134118301172172, 
-.192123665432045178897238849952, -.540265562572632296551266978691,
 .143620150765492549865289851418, 5.76672914073797259424119406295, 
-0.633944990402153813757581223306e-2, 0.626676247850351395534029175476e-2, 
.180009326296545764367708277240, -0.948059313982808533705260942973e-1, 
.175875688555801634799614555557, -0.633944990402153813757581223306e-2, 
.521502845604103687050998268457, -.969933707199744163847011746855, 
-.175103699996695180664621673847, 0.874409142585676611783165499024e-1,
 -.170820256565046134118301172172, 0.626676247850351395534029175476e-2,
 -.969933707199744163847011746855, 0.443320454485617625090646778007e-3), 6, 6)


cc <- chkerrors(loglik, param, y, x, truese)
print(cc)