# KOMMANDOER TIL FORELESNINGENE FOR UKE 14 

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


# EKSEMEL 1 FRA HEFTET TIL STORVIK: NEDBRYTING AV MYONER

# Leser inn dataene:

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)

# Histogram av dataene:
hist(x)

# 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 for eksempel 1
# (jf. sidene 10-14 i heftet til Storvik)

# Bruker startverdi alpha0=-0.5 og plotter kvadratisk 
# tiln?rming til log-likelihooden  
# [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 (svarer til kommandoene p? side 13
# i heftet til Storvik)

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)
}

# Estimat med standardfeil (basert p? observert informasjon):
cat("Estimat:", a, "Standardfeil:",sqrt(-1/d2loglik),"\n")



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

# EKSEMEL 2 FRA HEFTET TIL STORVIK: REGN I ILLINOIS


# Leser inn dataene:
x=
scan("http://www.uio.no/studier/emner/matnat/math/STK2120/v16/illrain.txt",na.strings="*")
x=x[!is.na(x)]

# lager et histogram:
hist(x)

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

# Beregner log-likelihooden og plotter den 
# (jf. ekstraoppgave 7)
# (svarer til kommandoene p? side 44 i heftet til Storvik)
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]))
persp(alpha,lambda,loglik,theta=330,phi=45,shade=1,
      zlab="loglik")
contour(alpha,lambda,loglik)

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

# Vi lager en 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) og 
# estimert kovariansmatrise COV 
#(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)


# Estimat med standardfeil:
cat("Estimat alpha:", theta[1], "Standardfeil:",
   sqrt(diag(COV))[1],"\n")
cat("Estimat lambda:", theta[2], "Standardfeil:",
   sqrt(diag(COV))[2],"\n")