# Kommandoer til ekstraoppgave 10 i Storviks hefte

# Leser inn dataene 
#(x="heart rate", y="oxygen uptake", begge som prosent av maksimal verdi)
x=c(43.5,44.0,44.0,44.5,44.0,45.0,48.0,49.0,49.5,51.0,54.5,57.5,57.7,61.0,63.0,72.0)
y=c(22.0,21.0,22.0,21.5,25.5,24.5,30.0,28.0,32.0,29.0,38.5,30.5,57.0,40.0,58.0,72.0)
n=length(x)


# Punkt a):

# Vanlig lineær regresjon
fit=lm(y~x)
summary(fit)

# Sjekk av normalitet
par(mfrow=c(1,2))
hist(fit$residuals)
qqnorm(fit$residuals)
par(mfrow=c(1,1))



# Punkt b)
# Bootstapping for regresjon ved å resample residualene

beta0.hat=fit$coef[1]
beta1.hat=fit$coef[2]
epsilon=fit$residuals
beta0.star=rep(NA,B)
beta1.star=rep(NA,B)
B=1000
for (b in 1:B)
{
x.star=x
epsilon.star=sample(epsilon,n,replace=T)
y.star=beta0.hat+beta1.hat*x.star+epsilon.star
fit.star=lm(y.star~x.star)
beta0.star[b]=fit.star$coef[1]
beta1.star[b]=fit.star$coef[2]
}

# Estimert skjevhet og standardfeil:
bias.beta0=mean(beta0.star)-beta0.hat
se.beta0=sd(beta0.star)
bias.beta1=mean(beta1.star)-beta1.hat
se.beta1=sd(beta1.star)


# Punkt c)
# Bootstapping for regresjon ved å resample (xi,yi)-ene
beta0.star=rep(NA,B)
beta1.star=rep(NA,B)
B=1000
for (b in 1:B)
{
ind=sample(1:n,n,replace=T)
x.star=x[ind]
y.star=y[ind]
fit.star=lm(y.star~x.star)
beta0.star[b]=fit.star$coef[1]
beta1.star[b]=fit.star$coef[2]
}

# Estimert skjevhet og standardfeil:
bias.beta0=mean(beta0.star)-beta0.hat
se.beta0=sd(beta0.star)
bias.beta1=mean(beta1.star)-beta1.hat
se.beta1=sd(beta1.star)