### Failure times of an engine component with different levels of corrosion.
### MLE's of a and b in the model

### ref: An Introduction to Statistical Modeling of Extreme Values
### Stuart Coles, pp. 38-43

library(ismev)
library(stats4)

# data from the ismev library

data(engine)

ft = engine$Time    # failure time
n = length(ft)      # sample size
w = engine$Corrosion   # corrosion measurement

plot(w,ft)

# MLE

# solve the nonlinear equation 

# minus the log likelihood

ll = function(a,b) { 
  -n*log(a) - b*sum(log(w)) + a*sum((w^b)*ft)
}

model.mle = mle(minuslog=ll,start=list(a=1,b=1)) 
model.mle

# plot fitted model

a.mle = coef(model.mle)[1]
a.mle
b.mle = coef(model.mle)[2]
b.mle

w.index = seq(min(w),max(w),0.01)

ft.fit = (a.mle^-1)*w.index^(-b.mle)

plot(w,ft,xlab="Corrosion Level",ylab="Lifetime",main="Fitted Model")
lines(w.index,ft.fit,type="l",col=3)

# value of the log-likelihood

ll.value = -ll(a.mle,b.mle)
ll.value

# CIs for a and b

a.mle

I.O = matrix(c( n*a.mle^-2, sum(w^b.mle*ft*log(w)),
sum(w^b.mle*ft*log(w)), a.mle*sum(w^b.mle*ft*log(w)^2) ), c(2,2))
I.O 

I.O.inv = solve(I.O)  # produces the matrix inverse
I.O.inv

a.mle.se = sqrt(I.O.inv[1,1])
a.mle.se

conf.level = 0.95
cv = qnorm(1-(1-conf.level)/2)

a.ci = c(a.mle - cv*a.mle.se,a.mle + cv*a.mle.se)
a.ci

b.mle

b.mle.se = sqrt(I.O.inv[2,2])
b.mle.se

b.ci = c(b.mle - cv*b.mle.se,b.mle + cv*b.mle.se)
b.ci

# Confidence Interval for phi = E[T|W=w0] = a^-1*w0^-b at some fixed value of w = w0
# using the delta method

w0 = 2

phi.mle = a.mle^-1*w0^(-b.mle)
phi.mle

delta.phi = matrix(c(-a.mle^-2*w0^-b.mle, -a.mle^-1*w0^b.mle*log(w0)),c(2,1))
delta.phi

delta.phi.T = t(delta.phi)   # transpose
delta.phi.T

V = I.O.inv

var.phi = delta.phi.T %*% V %*% delta.phi
var.phi

se.phi = sqrt(var.phi)
se.phi

phi.ci = c(phi.mle - cv*se.phi, phi.mle + cv*se.phi)
phi.ci

# Test b = 0, i.e., T ~ Exp(a)

# Fit the exponential model

a.mle.small = n/sum(ft)
a.mle.small

# plot the model

w.index = seq(min(w),max(w),0.01)

ft.fit.small = a.mle.small^-1

plot(w,ft,xlab="Corrosion Level",ylab="Lifetime",main="Fitted Models")
lines(w.index,ft.fit,type="l",col=3)
abline(h=ft.fit.small,col=4)

# value of the log-likelihood

ll.value.small = n*log(a.mle.small) - a.mle.small*sum(ft)
ll.value.small

# Deviance: D =  -2*log(Phi) = 2(ll.value - ll.value.small)

D = 2*(ll.value - ll.value.small)
D

# compare the Chi-Square 0.95 df = 1 value

cv = qchisq(0.95, df=1)
cv

# since the D > cv, we reject H0: b = 0