# KOMMANDOER TIL FORELESNINGENE FOR UKE 15 

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


# EKSEMEL 2 FRA HEFTET TIL STORVIK: REGN I ILLINOIS
#(se ogsaa ekstraoppgave 7)

# Vi vil bruke likelihood ratio testen til aa teste 
# nullhypotesen (H0) at alpha=1 (dvs. at datene er eksponensial fordelt)
# mot den alternative hypotesen (Ha) at alpha er forskjellig fra 1 

# Leser inn dataene:

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


# Vi lager en liten lkke som beregner ML-estimatene ved 
# Newton-Raphsons metode og tilhrende maksimumsverdi av log-likelihooden
# (jf. R-kommandoer til forelesningene for uke 14)

n = length(x)
sumx=sum(x)
sumlogx = sum(log(x));
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))
}

# Vi merker oss at "loglik" er maksimumsverdien av log-likelihooden

# Vi beregner ogs ML-estimatet for lambda under H0 med tilhoerende 
# maksimumsverdien av log-likelihooden:
lambda0=n/sumx
loglik0=n*log(lambda0)-lambda0*sumx
# Endelig beregner vi -2*log(sannsynlighetskvoten):

2*(loglik-loglik0)

# Vi finner at -2*log(sannsynlighetskvoten)=146.3. 
# Denne verdien skal sammenlignes med kji-kvadratfordelingen 
# med 2-1=1 frihetsgrad. Det gir en P-verdi tilnaermet lik 0, 
# saa H0 forkastes