# KOMMANDOER TIL OM BOOTSTRAP TIL FORELESNINGENE FOR UKE 19 

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


# PARAMETRISK BOOTSTRAP: SKJEVHET OG STANDARDFEIL
# EKSEMPEL: REGN I ILLINOIS 
# SE EKSEMPEL 2 FRA OPTIMERINGSHEFTET TIL STORVIK 
# SE OGSAA 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 foerst moment estimatene (brukes som startverdier for ML)
lambda.hat=mean(x)/var(x)
alpha.hat=mean(x)^2/var(x)

# Vi vil studere skjevhet og standardfeil for 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" med 
# momentestimatene som startverdi
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]
}

# Finner skjevhet og standardfeil for estimatene:
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 14 standardfeilene
# 0.0338 og 0.247 for alpha.ml og lambda.ml ved  aa bruke den
# tilnaermete fordelingen til ML-estimatene.

# Vi ser paa histogrammene for bootstrap-estimatene 
# (som gir en tilnaerming til fordelingen estimatene)
par(mfcol=c(2,1))
hist(alpha.ml.star,20,main="ML alpha")
hist(lambda.ml.star,20,main="ML lambda")
par(mfcol=c(1,1))




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

# FORVENTNING OG MEDIAN VED IKKE-PARAMETRISK BOOTSTRAP
# EKSEMPEL: VEKT?KNING I LOEPET AV ET SVANGERSKAP

# Vi ser paa dataene fra oppgave 8.51 i D&B. Her er det gitt 
# vekt?kning (i pund) i loepet av svangerskapet for 68 kvinner

# Vi leser inn dataene:
x=c(25,14,20,38,21,22,36,38,35,37,35,24,31,28,25,32,23,30,39,26, 38,20,21,11,35,42,31,25,59,23,43,38,21,76,22,26,10,19,25,25,15, 31,34,36,35,33,24,44,35,43,7,32,25,27,31,14,25,16,25,47,35,-14, 65,40,35,45,27,24)

n=length(x)

# Histogram av dataene:
hist(x)

# Vi er interessert i  estimere forventet vektoekning og median vektoekning

# Gjennomsnittet er
forv.hat=mean(x)

# De empiriske medianen er
med.hat=median(x)

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

# Bootstrap estimat for skjevhet
bias.forv=mean(forv.star)-forv.hat
bias.med=mean(med.star)-med.hat

# Bootstrap estimat for standardfeil:
se.forv=sd(forv.star)
se.med=sd(med.star)

# Vi ser paa histogrammene for bootstrap-estimatene 
# (som gir en tilnaerming til fordelingen estimatene)
par(mfcol=c(2,1))
hist(forv.star,20,main="Forventning")
hist(med.star,20,main="Median")
par(mfcol=c(1,1))

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

# BIVARIATE DATA (se avsnitt 5 notatet til Storvik om bootstrap;
# vi noeyer oss med aa se paa ikke-parametrisk bootstrap)

# Leser inn data om kastanjetraer
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)




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

# 95% KONFIDENSINTERVALL FOR MEDIAN VED IKKE-PARAMETRISK BOOTSTRAP
# EKSEMPEL: VEKT?KNING I LOEPET AV ET SVANGERSKAP

# Vi ser videre paa dataene fra oppgave 8.51 i D&B om
# vekt?kning (i pund) i loepet av svangerskapet for 68 kvinner

# Vi ser her paa medianen (forventnigen er tilsvarende)

# Vi ser paa dataene fra oppgave 8.51 i D&B. Her er det gitt 
# vekt?kning (i pund) i loepet av svangerskapet for 68 kvinner

# Vi gjentar noen av kommandoene slik at eksempelet blir fullstendig

# Vi leser inn dataene:
x=c(25,14,20,38,21,22,36,38,35,37,35,24,31,28,25,32,23,30,39,26, 38,20,21,11,35,42,31,25,59,23,43,38,21,76,22,26,10,19,25,25,15, 31,34,36,35,33,24,44,35,43,7,32,25,27,31,14,25,16,25,47,35,-14, 65,40,35,45,27,24)

n=length(x)

# De empiriske medianen er
med.hat=median(x)

# Vi bruker 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)
}

# Vi er interessert i  bestemme konfidensintervall for medianen
# Det er mange mulige konfidensintervaller

# Ser foerst paa konfidensintervall basert paa tilnaermet normalitet 
se.med=sd(med.star)
ki.med.norm=c(med.hat-1.96*se.med,med.hat+1.96*se.med)

# Ser saa paa standard ("basic") bootstrap konfidensintervall 
# (som er det som diskuteres i avsnitt 4 i Storviks hefte)

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

# Vi ser endelig paa percentil intervallet
# (som er det som diskuteres i avsnitt 8.5 i D&B)
ki.med.perc=c(med.sort[k1],med.sort[k2])
 


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

# BRUK AV BIBILOTEKET boot

# Vi har saa langt programert bootstrap (siden det gir forstaaelse 
# i hvordan metoden virker). Men beregningene blir enklere hvis vi 
# bruker bibiloteket boot.
# Vi viser framgangsmten for eksempelet med vektoekning 

library(boot)

# Vi gjentar noen av kommandoene slik at eksempelet blir fullstendig

# Vi leser inn dataene:
x=c(25,14,20,38,21,22,36,38,35,37,35,24,31,28,25,32,23,30,39,26, 38,20,21,11,35,42,31,25,59,23,43,38,21,76,22,26,10,19,25,25,15, 31,34,36,35,33,24,44,35,43,7,32,25,27,31,14,25,16,25,47,35,-14, 65,40,35,45,27,24)

n=length(x)

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

# Saa bruker vi funksjonen boot og ser paa resultatet:
boot.res=boot(x,median.boot,1000)
print(boot.res)
plot(boot.res)

# Saa finner vi ulike bootstrap konfidensintervall
boot.ci(boot.res,conf=0.95,type=c("norm","basic","perc"))

# Konfidensintervallet "norm" er her korrigert for skjevheten,
# saa i eksemplet ovenfor svaer dette til foelgende: 
ki.med.norm-(mean(med.star)-med.hat)