# KOMMANDOER TIL FORELESNINGENE FOR UKE 11 

# ========================================


# EKSEMEL 1 FRA HEFTET TIL STORVIK: NEDBRYTING AV MYONER

# Leser inn dataene (jf. side 1 i heftet til Storvik)

x=c(0.41040018,0.91061564,-0.61106896,0.39736684,0.37997637,0.34565436 ,
0.01906680,-0.28765977,-0.33169289,0.99989810,-0.35203164,0.10360470 ,
0.30573300,0.75283842,-0.33736278,-0.91455101,-0.76222116,0.27150040 ,
-0.01257456,0.68492778,-0.72343908,0.45530570,0.86249107,0.52578673 ,
0.14145264,0.76645754,-0.65536275,0.12497668,0.74971197,0.53839119)

# --------

# Beregner log-likelihooden og plotter den 
#(jf. sidene 1-2 i heftet til Storvik)

n=length(x)
alpha=seq(-1,1,0.01)
loglik=rep(0,length(alpha))
for (i in 1:length(alpha))
loglik[i]=sum(log(1+alpha[i]*x))-n*log(2)
plot(alpha,loglik,type='l')

#----------

# Illustrasjon av Newton-Raphsons metode

# Bruker startverdi alpha0=-0.5 og plotter kvadratisk tiln?rming 
# til log-likelihooden med [jf. formel (4.1) i heftet til Storvik] 

alpha0=-0.5
loglik0=sum(log(1+alpha0*x))-n*log(2)
dloglik0=sum(x/(1+alpha0*x))
d2loglik0=-sum((x/(1+alpha0*x))^2)
kvloglik0=loglik0+dloglik0*(alpha-alpha0)+0.5*d2loglik0*(alpha-alpha0)^2
lines(alpha,kvloglik0,col='red')

# Den maksimerende verdien for den kvadratiske tiln?rmingen er 
# alpha1=0.085 (se nedenfor). Vi plotter ny kvadratisk tiln?rming 
# til log-likelihooden ut fra denne verdien

alpha1=0.085
loglik1=sum(log(1+alpha1*x))-n*log(2)
dloglik1=sum(x/(1+alpha1*x))
d2loglik1=-sum((x/(1+alpha1*x))^2)
kvloglik1=loglik1+dloglik1*(alpha-alpha1)+0.5*d2loglik1*(alpha-alpha1)^2
lines(alpha,kvloglik1,col='blue')

# Den maksimerende verdien for den nye kvadratiske tiln?rming er 
# alpha2=0.516 (se nedenfor). Vi lager enda en ny kvadratisk tiln?rming 
# ut fra denne, og finner ny maksimerende verdi. Vi fortsetter slik til
# alpha-verdiene ikke endrer seg mer enn epsilon (dvs et gitt krav til
# n?yaktighet. Dette er Newton-Raphsons metode.

#-----------------


# Vi lager en liten l?kke som gir ML-estimatet ved Newton-Raphsons metode 

eps=0.000001
diff=1
a=-0.50
i=0
loglik=sum(log(1+a*x))-n*log(2)
dloglik=sum(x/(1+a*x))
d2loglik=-sum((x/(1+a*x))^2)
cat("iterasjon=",i,"alpha=",a,"loglik=",loglik,"dloglik=",dloglik,"d2loglik=", d2loglik,"\n")
while(diff>eps)
{
i=i+1
a.old=a
a=a-dloglik/d2loglik
loglik=sum(log(1+a*x))-n*log(2)
dloglik=sum(x/(1+a*x))
d2loglik=-sum((x/(1+a*x))^2)
cat("iterasjon=",i,"alpha=",a,"loglik=",loglik,"dloglik=",dloglik,"d2loglik=", d2loglik,"\n")
diff=abs(a-a.old)
}

# =======================

# EKSEMEL 2 FRA HEFTET TIL STORVIK: REGN I ILLINOIS
#(se ogs? ekstraoppgave 6)

# Leser inn dataene:

x=
scan("http://www.uio.no/studier/emner/matnat/math/STK2120/v13/illrain.dat",
na.strings="*")
x=x[!is.na(x)]

#------------------

# Beregner log-likelihooden og plotter den (jf. ekstraoppgave 6)

alpha = seq(0.35,0.55,0.005)
lambda = seq(1,3,0.05)
loglik = matrix(nrow=length(alpha),ncol=length(lambda))
n = length(x)
sumx=sum(x)
sumlogx = sum(log(x));
for(i in 1:length(alpha))
for(j in 1:length(lambda))
loglik[i,j] = n*alpha[i]*log(lambda[j])+(alpha[i]-1)*sumlogx-
lambda[j]*sumx-n*log(gamma(alpha[i]))
par(mfrow=c(1,2))
persp(alpha,lambda,loglik,theta=330,phi=45,shade=1,zlab="loglik")
contour(alpha,lambda,loglik)

#------------------------

# Vi lager en liten l?kke som gir ML-estimatet ved Newton-Raphsons metode 
# Vi bruker momentestimatene som startverdier. Merk at digamma-funksjonen 
# gir den deriverte av logaritmen til gamma-funksjonen, mens 
# trigamma-funksjonen gir den andre deriverte 

eps=0.000001
diff = 1
alpha = mean(x)^2/var(x)
lambda=mean(x)/var(x)
theta = c(alpha,lambda)
i=0
loglik=n*alpha*log(lambda)+(alpha-1)*sumlogx-
lambda*sumx-n*log(gamma(alpha))
cat("iterasjon=",i,"alpha=",alpha,"lambda=",lambda,"loglik=",loglik,"\n")

while(diff>eps)
{
i=i+1
theta.old=theta
s=c(n*log(lambda)+sumlogx-n*digamma(alpha),n*alpha/lambda-sumx)
J=matrix(c(n*trigamma(alpha),-n/lambda,-n/lambda,n*alpha/lambda^2),2,2)
theta=theta+solve(J)%*%s
alpha=theta[1]
lambda=theta[2]
loglik=n*alpha*log(lambda)+(alpha-1)*sumlogx-lambda*sumx-n*log(gamma(alpha))
cat("iterasjon=",i,"alpha=",alpha,"lambda=",lambda,"loglik=",loglik,"\n")
diff=sum(abs(theta-theta.old))
}

#------------------------

# Beregner observert informasjonsmatrise J (som her er lik Fishers
# informasjonsmatrise), estimert kovariansmatrise COV og standardfeilene 
# til estimatene (jf. sidene 9-10 og 15 i heftet til Storvik) 

J=matrix(c(n*trigamma(alpha),-n/lambda,-n/lambda,n*alpha/lambda^2),2,2)
COV=solve(J)
sqrt(diag(COV))