# KOMMANDOER TIL FORELESNINGENE FOR UKE 19 


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

# EKSEMPEL: REGN I ILLINOIS (MOMENTESTIMATER OG KONFIDENSINTERVALL)
# SE EKSEMPEL 2 FRA OPTIMERINGSHEFTET TIL STORVIK 
# SE OGS? SIDENE 6-7 I STORVIKS HEFTE OM BOOTSTRAP

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

# Vi finner moment estimatene:
lambda.hat=mean(x)/var(x)
alpha.hat=mean(x)^2/var(x)

# Vi vil f?rst studere skjevhet og standardfeil for momentestimatene
# Kommandoene blir tilsvarende som de p? side 7 i bootstrap heftet.
B=1000
lambda.star=rep(NA,B)
alpha.star=rep(NA,B)
for (b in 1:B)
{
x.star=rgamma=rgamma(n,shape=alpha.hat,rate=lambda.hat)
lambda.star[b]=mean(x.star)/var(x.star)
alpha.star[b]=mean(x.star)^2/var(x.star)
}
bias.lambda=mean(lambda.star)-lambda.hat
se.lambda=sd(lambda.star)
bias.alpha=mean(alpha.star)-alpha.hat
se.alpha=sd(alpha.star)


# Vi vil s? se hvordan vi kan lage konfidensintervall for med 
# utgangspunkt i momentestimatene. Til illustrasjon ser vi p? 
# konfidensintervall for lambda (det blir tilsvarende for alpha)

# En mulighet er ? ta utgangspunkt i at momentestimatene er tiln?rmet 
# normalfordelt. Det kan vi sjekke ved ? se p? fordelingen til 
# lambda.star(som gir en tiln?rmelse til fordelingen til lambda.hat):
par(mfrow=c(1,2))
hist(lambda.star,20)
qqnorm(lambda.star)
par(mfrow=c(1,1))


# Et 95% konfidensintervall basert p? normalitet f?s ved:
ki.lambda.normal=c(lambda.hat-1.96*se.lambda, lambda.hat+1.96*se.lambda)

# En annen mulighet er ? beregne et standard bootstrap konfdensintervall
# slik det er beskrevet p? side 8 i heftet til Storvik. 
# (Dette intervallet forutsetter ikke tiln?rmet normalitet.)

# Histogram for lambda.star-lambda.hat:
hist(lambda.star-lambda.hat,20)

# Beregner empiriske 2.5% og 97.5% kvantiler
k1=as.integer(0.025*B)
k2=as.integer(0.975*B)
lambda.sort=sort(lambda.star)
deltaL=lambda.sort[k1]-lambda.hat
deltaU=lambda.sort[k2]-lambda.hat

# Standard bootstrap konfidensintervall:
ki.lambda.boot=c(lambda.hat-deltaU,lambda.hat-deltaL)




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

# EKSEMPEL: REGN I ILLINOIS (SAMMENLIGNING AV MOMENT- OG ML-ESTIMATER)

# Vi vil s? se p? ML-estimatene
# 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":
par=optim(c(alpha.hat,lambda.hat),negloglik,data=x)$par
alpha.ml=par[1]
lambda.ml=par[2]


# Bootstrap:
B=1000
lambda.ml.star=rep(NA,B)
alpha.ml.star=rep(NA,B)
for (b in 1:B)
{
x.star=rgamma(n,shape=alpha.ml,rate=lambda.ml)
par=optim(c(alpha.ml,lambda.ml),negloglik,data=x.star)$par
alpha.ml.star[b]=par[1]
lambda.ml.star[b]=par[2]
}
bias.lambda.ml=mean(lambda.ml.star)-lambda.ml
se.lambda.ml=sd(lambda.ml.star)
bias.alpha.ml=mean(alpha.ml.star)-alpha.ml
se.alpha.ml=sd(alpha.ml.star)


# Til sammenligning fant vi i uke 11 fant vi standardfeilene
# 0.0338 og 0.247 for alpha.ml og lambda.ml ved ? bruke den
# tiln?rmete fordelingen til ML-estimatene.

# Vi vi sammenligne fordelingen for momentestimatene og ML-estimatene.
# Vi ser da p? histogrammene for bootstrap-estimatene 
# (som gir en tiln?rming til fordelingen estimatene)
par(mfcol=c(2,2))
hist(alpha.star,20,xlim=c(min(alpha.star),max(alpha.star)), 
     main="Moment alpha")
hist(alpha.ml.star,20,xlim=c(min(alpha.star),max(alpha.star)),
     main="ML alpha")
hist(lambda.star,20,xlim=c(min(lambda.star),max(lambda.star)),
     main="Moment lambda")
hist(lambda.ml.star,20,xlim=c(min(lambda.star),max(lambda.star)),
     main="ML lambda")
par(mfcol=c(1,1))



# Et 95% konfidensintervall for lambda basert p? normalitet f?s ved:
ki.lambda.ml.normal=c(lambda.ml-1.96*se.lambda.ml, 
                       lambda.ml+1.96*se.lambda.ml)

# En standard bootstrap konfidensintervall f?s ved:
lambda.ml.sort=sort(lambda.ml.star)
deltaL.ml=lambda.ml.sort[k1]-lambda.ml
deltaU.ml=lambda.ml.sort[k2]-lambda.ml

# Standard bootstrap konfidensintervall:
ki.lambda.ml.boot=c(lambda.ml-deltaU.ml,lambda.ml-deltaL.ml)



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

# 95% KONFIDENSINTERVALL FOR MEDIAN VED IKKE-PARAMETRISK BOOTSTRAP

# Vi ser p? dataene for behandlingsgruppen i eksempel 1 i Storviks notat. 
# Der er det gitt levetider for 7 behandlede mus.

# Vi leser inn dataene:
x=c(94,197,16,38,99,141,23)
n=length(x)

# Vi er interessert i ? bestemme konfidensintervall for median levetid
# De empiriske medianen er
med.hat=median(x)

# Vi vil bruke ikke-parametrisk bootstrap 
B=1000
med.star=rep(NA,B) 
for (b in 1:B)
{
x.star=sample(x,n,replace=T)
med.star[b]=median(x.star)
}

# Histogram for med.star-med.hat:
hist(med.star-med.hat)

# Beregner empiriske 2.5% og 97.5% kvantiler
k1=as.integer(0.025*B)
k2=as.integer(0.975*B)
med.sort=sort(med.star)
deltaL=med.sort[k1]-med.hat
deltaU=med.sort[k2]-med.hat

# Standard bootstrap konfidensintervall:
ki.med.basic=c(med.hat-deltaU,med.hat-deltaL)


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

# BIVARIATE DATA (se avsnitt 5 notatet til Storvik om bootstrap;
# vi n?yer oss med ? se p? ikke-parametrisk bootstrap)

# Leser inn data om kastanjetr?r
kastanje=read.table("http://www.uio.no/studier/emner/matnat/math/STK2120/v13/chestnut.txt", header=T)
n=dim(kastanje)[1]

# Plotter data og beregner korrelasjonen:
plot(kastanje$Age,kastanje$DBH)
rho.hat=cor(kastanje$Age,kastanje$DBH)

# Bootstrap:
B=1000
rho.star=rep(NA,B)
for (b in 1:B)
{
ind=sample(1:n,n,replace=T)
age.star=kastanje$Age[ind]
dbh.star=kastanje$DBH[ind]
rho.star[b]=cor(age.star,dbh.star)
}
bias.rho=mean(rho.star)-rho.hat
se.rho=sd(rho.star)

# Histogram for rho.star:
hist(rho.star,20)

# Beregner empiriske 2.5% og 97.5% kvantiler
k1=as.integer(0.025*B)
k2=as.integer(0.975*B)
rho.sort=sort(rho.star)
deltaL=rho.sort[k1]-rho.hat
deltaU=rho.sort[k2]-rho.hat

# Standard bootstrap konfidensintervall:
ki.rho.basic=c(rho.hat-deltaU,rho.hat-deltaL)


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

# BRUK AV BIBILOTEKET boot

# Vi har s? langt programert bootstrap (siden det gir forst?else 
# i hvordan metoden virker). Men beregningene blir enklere hvis vi 
# bruker bibiloteket boot.
# Vi viser framgangsm?ten for det siste eksempelet

library(boot)

# Vi definerer en funksjon som regner ut korrelasjonen for et utvalg:
cor.boot=function(x,ind){cor(x[ind,1],x[ind,2])}

# S? bruker vi funksjonen boot og ser p? resultatet
boot.res=boot(kastanje,cor.boot,1000)
print(boot.res)
plot(boot.res)

# S? finner vi standard bootstrap konfidensintervall
boot.ci(boot.res,conf=0.95,type='basic')