# KOMMANDOER TIL OM BOOTSTRAP TIL FORELESNINGENE FOR UKE 17 

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

# MONTE CARLO INTEGRASJON

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

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

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


# Monte Carlo integrasjon er s?rlig nyttig for integraler av h?yere 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 opph?yd i 10. Men metoden er generell og kan 
# brukes tilsvarende for andre funksjoner av u1;?,u10.


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

# ESTIMERING AV FORVENTNING (=MEDIAN) I NORMALFORDELINGEN

# Vi ser p? 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 estimerer f?rst 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 n?dvendig her, men gir 
# en metode som # kan brukes generelt.) 
# Vi trekker da B=1000 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 v?re et estimat for standardfeilen:

B=1000
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 p? kalles parametrisk bootstrap


# Vi vil s? 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 f?lge oppgave 7.48 i D&B kan en vise at variansen til den empiriske
# medianen er tiln?rmet lik 1/[4*n*f(my)^2]. Det gir standardfeilen
se2=sqrt(1/(4*n*dnorm(my2,my1,s)^2))

# Her vil bootstrap estimatet v?re n?yaktigere enn 
# det tiln?rmete resultatet

# Vi kan ogs? bruke bootstrap metoden til ? bestemme skjevheten 
# til den empiriske medianen:
bias.my2=mean(my2.star)-my2