# KOMMANDOER TIL FORELESNINGENE FOR UKE 14 

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


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

# Vi vil bruke sannsynlighetskvotetesten til å teste 
# nullhypotesen (H0) at alpha=1 (dvs. at datene er eksponensial fordelt)
# mot den alternative hypotesen at alpha er forskjellig fra 1 

# Leser inn dataene:

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


# Vi lager en liten løkke som beregner ML-estimatene ved 
# Newton-Raphsons metode og tilhørende maksimumsverdi av log-likelihooden
# (jf. R-kommandoer til forelesningene for uke 11)

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 tilhørende 
# 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 tilnærmet lik 0, 
# så H0 forkastes