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

    # -------------------------------------------
    # Nuclear power station example - Example 5.2
    # from Brazzale, Davison and Reid (2007)
    #
    # This file does the OLS fit
    # -------------------------------------------
library(marg)
data(nuclear)

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



#===============================
 loglik <- function( beta, y, x)
#===============================
{
    # -------------------------------------------------
    # Log-likelihood function for nuclear power station 
    # data - model with t_4 errors
    # -------------------------------------------------

    mu = beta[1] + beta[2]*x[,1] + beta[3]*x[,2] + 
         beta[4]*x[,3] + beta[5]*x[,4] + beta[6]*x[,5] + 
         beta[7]*x[,6]
    z = (y-mu)/(exp(beta[8]))
    nu = 4
    const = gamma((nu+1)/2) / ( sqrt(pi*nu) * gamma(nu/2) )
    sum(  -log(const) + (nu+1)/2 * log(1+z*z/nu)  + (beta[8]))
}



inv4 = matrix( c(
13.4776796619498652955188148026, -.173706259102035477656680169708, 
-.252566829158921953518192158193, -.114436976716955879619098751902,
 -0.368678054796215373145712644146e-1, .112077606879754435719227309473, 
-.290601515480106078142662837366, -0.474211377016611442823343670749e-1, 
-.173706259102035477656680169633, 0.231401068130629530225019627986e-2, 
0.247647221515522905252009335222e-2, 0.160888742182657422048241796684e-2,
 0.530946592054595024440720341280e-3, -0.145915640795081244978694093232e-2,
 0.392348789314640063755988887045e-2, 0.425033423118891061744612648788e-3, 
-.252566829158921953518192158992, 0.247647221515522905252009336406e-2,
 0.128754364522310957891375982543e-1, 0.823388878068581332576259508579e-3,
 -0.587439210264119018569704316766e-4, -0.215757000798419440777080690620e-2,
 0.356078304891461100389200827295e-2, 0.281577543257272810810586143237e-2, 
-.114436976716955879619098751880, 0.160888742182657422048241796731e-2,
 0.823388878068581332576259499725e-3, 0.600059571638512657802326348613e-2,
 -0.485666306710116338651254689734e-3, -0.177477683078095798327324060226e-2,
 0.425774410065529368023684371650e-2, -0.484778221246383593856309018541e-3, 
-0.368678054796215373145712642485e-1, 0.530946592054595024440720339293e-3, 
-0.587439210264119018569704366369e-4, -0.485666306710116338651254691273e-3,
 0.287954432722052246711129687359e-2, -0.170760051877283621237095019392e-3,
 0.133176353787282808350421292008e-2, -0.456668308278013962114081572818e-3, 
.112077606879754435719227309682, -0.145915640795081244978694093580e-2,
 -0.215757000798419440777080690157e-2, -0.177477683078095798327324060453e-2, 
-0.170760051877283621237095021484e-3, 0.183214014769894465030982069447e-2,
 -0.306958406719972499185792652931e-2, -0.483685350219924762973513713256e-3, 
-.290601515480106078142662837227, 0.392348789314640063755988887033e-2, 
0.356078304891461100389200825170e-2, 0.425774410065529368023684371570e-2,
 0.133176353787282808350421292349e-2, -0.306958406719972499185792652324e-2,
 0.101027824732158330836886675905e-1, 0.897837574522931830220730459133e-3, 
-0.474211377016611442823343669619e-1, 0.425033423118891061744612647594e-3, 
0.281577543257272810810586142741e-2, -0.484778221246383593856309019000e-3, 
-0.456668308278013962114081572642e-3, -0.483685350219924762973513712051e-3,
 0.897837574522931830220730457534e-3, 0.289176609448593411251900924024e-1), 8, 8)

inv20 = matrix(c(8.83005931940038533786022853634, -.116468822098501006991321501764,
 -.136041735051859427963943685691, -0.600001600165361801583794005080e-1, 
-0.357045267518127838048641079646e-1, 0.756474371561481002518785815941e-1,
 -.213142365349529496800533339384, -0.147535200156730565737130901852e-1, 
-.116468822098501006991321501985, 0.162366096638860605670888770285e-2, 
0.894070550895747884468000676932e-3, 0.823493429679244231679627534609e-3, 
0.495956420391905244157346810241e-3, -0.100869386895026876457161541304e-2, 
0.297445267685202460951073359249e-2, 0.172313187822842296019078247625e-3,
 -.136041735051859427963943683401, 0.894070550895747884468000643383e-3,
 0.114491457931507579788751452292e-1, 0.572693774620688527858306341745e-3,
 0.785151023876129717768274649464e-4, -0.124618630008801817487820454169e-2,
 0.151515989105656121633554082700e-2, 0.467293676297108858447608759992e-3, 
-0.600001600165361801583794006170e-1, 0.823493429679244231679627534727e-3, 
0.572693774620688527858306357035e-3, 0.451663727917507871219908159773e-2, 
-0.146935540287480977960541021233e-3, -0.107532982693974712108392165095e-2, 
0.295275505214610135268929763774e-2, 0.106194762764783133410611029367e-3, 
-0.357045267518127838048641081179e-1, 0.495956420391905244157346811299e-3,
 0.785151023876129717768274772538e-4, -0.146935540287480977960541020757e-3,
 0.295022710155771309917110250711e-2, -0.159337355321123599079866883189e-3,
 0.155469954471845648760043096958e-2, 0.191361263522302813722243028197e-4, 
0.756474371561481002518785816072e-1, -0.100869386895026876457161541163e-2,
 -0.124618630008801817487820455827e-2, -0.107532982693974712108392165024e-2, 
-0.159337355321123599079866882090e-3, 0.147432614563503990522107370501e-2, 
-0.240460203997957392199828647187e-2, -0.225572597109980925694893590537e-3, 
-.213142365349529496800533339875, 0.297445267685202460951073359381e-2,
 0.151515989105656121633554088796e-2, 0.295275505214610135268929763817e-2,
 0.155469954471845648760043096809e-2, -0.240460203997957392199828647513e-2,
 0.100869086576966872269179884770e-1, 0.474580555520938475248171751756e-3, 
-0.147535200156730565737130900909e-1, 0.172313187822842296019078246047e-3, 
0.467293676297108858447608762049e-3, 0.106194762764783133410611028614e-3, 
0.191361263522302813722243021770e-4, -0.225572597109980925694893589793e-3,
 0.474580555520938475248171748772e-3, 0.179558234645844641497391599644e-1), 8, 8)

inv99 = matrix(c(7.92395485779676530983594450766, -.105501606439438003275632065435,
 -.111352291360017709839544035395, -0.542559213410502125900021547776e-1, 
-0.322347733842449588364784307216e-1, 0.676363370862171912536061750017e-1, 
-.196422463160563773235510175282, -0.283712450554677707070244974177e-2, 
-.105501606439438003275632065269, 0.149442212294994356527673100421e-2, 
0.561071585865136341795431974414e-3, 0.742014244063100367579552556717e-3,
 0.442017876475681181930291552762e-3, -0.914408756690396217002815983483e-3, 
0.277267274748348564125116004615e-2, 0.340778508856527395631727724277e-4, 
-.111352291360017709839544037167, 0.561071585865136341795432000418e-3, 
0.111048406406759327623325742817e-1, 0.526400225916584535999686235951e-3, 
0.117849539183185855642156311891e-3, -0.988840983144702604622131502815e-3, 
0.105321895992369746338981559855e-2, 0.792956655554635804394962016343e-4, 
-0.542559213410502125900021546566e-1, 0.742014244063100367579552556408e-3, 
0.526400225916584535999686220565e-3, 0.433875367042646025996406595291e-2, 
-0.234398967244039730419850601027e-4, -0.990594326669215847603163879327e-3, 
0.281408726136356931147985837255e-2, 0.248529065923803999144599135805e-4, 
-0.322347733842449588364784306164e-1, 0.442017876475681181930291552029e-3, 
0.117849539183185855642156303364e-3, -0.234398967244039730419850603955e-4, 
0.287616819832309823469756577323e-2, -0.177392584216884460528715152787e-3, 
0.151028114910147832884170688216e-2, 0.969280845953607293018648343535e-5, 
0.676363370862171912536061750630e-1, -0.914408756690396217002815986042e-3, 
-0.988840983144702604622131485318e-3, -0.990594326669215847603163880853e-3, 
-0.177392584216884460528715154044e-3, 0.137012219495738129985434165435e-2, 
-0.225521636334432972341177332017e-2, -0.441320697196880282571818876658e-4, 
-.196422463160563773235510174949, 0.277267274748348564125116004611e-2, 
0.105321895992369746338981554740e-2, 0.281408726136356931147985837305e-2, 
0.151028114910147832884170688353e-2, -0.225521636334432972341177331527e-2, 
0.100827019694216368158889384532e-1, 0.994220809909936463589956290828e-4, 
-0.283712450554677707070244974211e-2, 0.340778508856527395631727725430e-4, 
0.792956655554635804394962004857e-4, 0.248529065923803999144599136452e-4,
 0.969280845953607293018648349699e-5, -0.441320697196880282571818876580e-4,
 0.994220809909936463589956293159e-4, 0.160733441417715361413067207832e-1), 8, 8)




param4 = c(-112951236, 1907705, 6475260, 2418897, 1439194, -598352, -2815051, -22009605)
param4 = param4 / 10000000

inv = inv4
for (j in 1:7) {
    for (k in (j+1):8) {
        inv[j,k] = inv[j,k] / sqrt(inv[j,j] * inv[k,k])
        inv[k,j] = inv[j,k]
    }
}
for (j in 1:8) {
    inv[j,j] = sqrt(inv[j,j])
}



nuc.t4 <- rsm( log(cost) ~ date + log(cap) + NE + CT + log(N) + PT,
                data=nuclear, dispersion = NULL, family=student(4))
summary(nuc.t4)


param = param4
truese = inv

attach(nuclear)
y = log(cost)
x = cbind(date,log(cap),NE,CT,log(N),PT)
detach(nuclear)

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