# Kommandoer til ekstraoppgave 9 i Geir Storviks hefte

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



# Punkt a):

# Ser p? konfidensintervall for alfa og lambda basert p? ML-estimatene og 
# parametrisk bootstrap. (Det blir tilsvarende for intervallene basert 
# p? momentestimatene, se R-kommandoene til forelesningene til uke 12)

# Vi definerer en funksjon som beregner minus log-likelihood:
negloglik=function(par,data)
{
alpha=par[1];lambda=par[2];n=length(data)
loglik=n*alpha*log(lambda)+(alpha-1)*sum(log(data))-
       lambda*sum(data)-n*log(gamma(alpha))
-loglik
}

# Finner ML-estimatene ved ? bruke kommandoen "optim" (kan v?re bedre 
# enn nlm for optimering i to dimensjoner).
# Vi bruker momentestimatene som startverdier
lambda.mom=mean(x)/var(x)
alpha.mom=mean(x)^2/var(x)
par=optim(c(alpha.mom,lambda.mom),negloglik,data=x)$par
alpha.hat=par[1]
lambda.hat=par[2]


# Parametrisk bootstrap:
B=1000
lambda.star=rep(NA,B)
alpha.star=rep(NA,B)
for (b in 1:B)
{
x.star=rgamma(n,shape=alpha.hat,rate=lambda.hat)
par=optim(c(alpha.hat,lambda.hat),negloglik,data=x.star)$par
alpha.star[b]=par[1]
lambda.star[b]=par[2]
}

# Standard (basic) bootstrap konfidensintervall for alpha f?s ved:
k1=as.integer(0.025*B)
k2=as.integer(0.975*B)
alpha.sort=sort(alpha.star)
deltaL.alpha=alpha.sort[k1]-alpha.hat
deltaU.alpha=alpha.sort[k2]-alpha.hat
ki.alpha.boot=c(alpha.hat-deltaU.alpha,alpha.hat-deltaL.alpha)
  
# Standard (basic) bootstrap konfidensintervall for lambda f?s ved:
lambda.sort=sort(lambda.star)
deltaL.lambda=lambda.sort[k1]-lambda.hat
deltaU.lambda=lambda.sort[k2]-lambda.hat
ki.lambda.boot=c(lambda.hat-deltaU.lambda,lambda.hat-deltaL.lambda)


# Punkt b)

# Vi beregner spredningskoeffisienten:
C.hat=sd(x)/mean(x)
C_i_empirisk=sqrt((n-1)/n)*C.hat

# Vi bruker ikke-parametrisk bootstrap 
B=1000
C.star=rep(NA,B) 
for (b in 1:B)
{
x.star=sample(x,n,replace=T)
C.star[b]=sd(x.star)/mean(x.star)
}

# Estimert skjevhet og standardfeil:
bias.C=mean(C.star)-C_i_empirisk
se.C=sd(C.star)

# Standard bootstrap konfidensintervall:
C.sort=sort(C.star)
deltaL=C.sort[k1]-C.hat
deltaU=C.sort[k2]-C.hat
ki.C.basic=c(C.hat-deltaU,C.hat-deltaL)