# KOMMANDOER TIL OM BOOTSTRAP TIL FORELESNINGENE FOR UKE 18 

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

# MONTE CARLO INTEGRASJON

# Vi vil bestemme integralet av exp(u^4) over intervallet [0,1]

# Vi beregner foerst integralet numerisk:
fx=function(x) exp(x^4)
integrate(fx,0,1)

# Vi beregner saa integralet (tilnaermet) ved Monte Carlo integrasjon:
B=10000
u=runif(B,0,1)
x=exp(u^4)
mean(x)


# Monte Carlo integrasjon er saerlig nyttig for integraler av hoeyere orden
# Vi beregner det 10-dimensjonale integralet av exp(u1^4+u2^4+...+u10^4)
B=10000
int=rep(NA,B)
for (b in 1:B)
{
u=runif(10,0,1)
int[b]=exp(sum(u^4))
}
mean(int)

# I dette er det 10-dimensjonale integralet lik det univariate 
# integralet ovenfor opphyd i 10. Men metoden er generell og kan 
# brukes tilsvarende for andre funksjoner av u1,....,u10.


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

# ESTIMERING AV FORVENTNING (=MEDIAN) I NORMALFORDELINGEN

# Vi ser paa dataene i eksempel 8.12 på side 404-405 i D&B. 
# Der er det gitt fettinnholdet for n=10 tilfeldig valgte hot dogs.

# Vi leser inn dataene:
fat=c(25.2,21.3,22.8,17.0,29.8,21.0,25.5,16.0,20.9,19.5)
n=length(fat)

# Vi estimerer foerst forventningen (=medianen) my ved gjennomsnittet 
my1=mean(fat)

# Vi estimerer standardfeilen til my1 ved
s=sd(fat)
se1=s/sqrt(n)

# Vi kan alternativt finne standardfeilen (= standardavviket for 
# estimatoren) ved simulering. (Det er ikke noedvendig her, men gir 
# en metode som kan brukes generelt.)
 
# Vi trekker da B=10000 ganger n=10 observasjoner fra normalfordelingen
# med forventning og standardavvik lik de estimerte. For hvert utvalg 
# estimerer vi forventningen, og til slutt bestemmer vi det empiriske 
# standardavviket. Det vil vaere et estimat for standardfeilen:

B=10000
my1.star=rep(NA,B)
for (b in 1:B)
{
x.star=rnorm(n,my1,s)
my1.star[b]=mean(x.star)
}
se1.b=sd(my1.star)

# Denne metoden  finne standardfeilen paa kalles parametrisk bootstrap


# Vi vil saa estimere my ved den empiriske medianen
my2=median(fat)

# Her er det vanskeligere  finne et uttrykk for standardfeilen.
# Men vi kan bruke bootstrap metoden:
B=10000
my2.star=rep(NA,B)
for (b in 1:B)
{
x.star=rnorm(n,my1,s)
my2.star[b]=median(x.star)
}
se2.b=sd(my2.star)

# I foelge oppgave 7.48 i D&B kan en vise at variansen til den empiriske
# medianen er tilnrmet lik 1/[4*n*f(my)^2]. Det gir standardfeilen
se2=sqrt(1/(4*n*dnorm(my2,my1,s)^2))

# Her vil bootstrap estimatet vaere nyaktigere enn 
# det tilnaermete resultatet


# Vi vil endelig se hvordan vi kan studere egenspapene til 
# den empiriske variasjonskoeffisienten 
cv=sd(fat)/mean(fat)

# Vi bruker bootstrap metoden:
B=10000
cv.star=rep(NA,B)
for (b in 1:B)
{
x.star=rnorm(n,my1,s)
cv.star[b]=sd(x.star)/mean(x.star)
}

# Vi estimerer skjevhet og standardfeil:
bias.cv=mean(cv.star)-cv
se_cv.b=sd(cv.star)



# EKSEMPEL: REGN I ILLINOIS 
# 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 foerst moment estimatene:
lambda.hat=mean(x)/var(x)
alpha.hat=mean(x)^2/var(x)

# Vi vil studere skjevhet og standardfeil for momentestimatene
# Kommandoene er tilsvarende som paa 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 saa 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" 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]
}
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 vi sammenligne fordelingen for momentestimatene og ML-estimatene.
# Vi ser da paa histogrammene for bootstrap-estimatene 
# (som gir en tilnaerming 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))



# STANDARDFEIL FOR MEDIAN VED IKKE-PARAMETRISK BOOTSTRAP

# Vi ser igjen paa dataene i eksempel 8.12 i D&B. 
# Der er det gitt fettinnholdet
# for n=10 tilfeldig valgte hot dogs.

# Vi leser inn dataene:
fat=c(25.2,21.3,22.8,17.0,29.8,21.0,25.5,16.0,20.9,19.5)
n=length(fat)

# Vi er interessert i  estimere medianen for fettinnholdet
# De empiriske medianen er
med=median(fat)

# Vi saa ovenfor paa parametrisk bootstrap
# (naar vi antok normalfordeling)


# Vi bruker naa ikke-parametrisk bootstrap. 
# Den eneste forskjellen er at vi da trekker x*-ene fra 
# den empiriske kumulative fordelingen, og det 
# svarer til trekke fra de observerte x-ene med tilbakelegging

med.i.star=rep(NA,B)
for (b in 1:B)
{
x.star=sample(fat,n,replace=T)
med.i.star[b]=median(x.star)
}
se.i.b=sd(med.i.star)