### bootstrap V & R 5.6

### variance of the mean of the normal distribution 
#  1000 samples from N(0,1), each of size 30

win.graph()
m <- 1000
res <- numeric(m)
for (i in 1:m) res[i] <- mean(rnorm(30))

print("standard deviation of 1000 means of samples of size 30")
print(sqrt(var(res)))                          # 1000*30 observations

readline()

# normal theory gives s/sqrt(n)
print("analytical result given by normal theory")
print(1/sqrt(30))

readline()
## distribution of means from repeated samples from N(0,1)


print("distribution of means")
hist(res, freq=FALSE)
lines(density(res))

readline()

# c.i.
print("confidence limits")
print(quantile(res, c(0.025, 0.975)))

readline()

##
## 300 samples from a sample of size 30.
##

sample30 <- rnorm(30)
print("sample of size 30")
print(sample30)

readline()

# how sample works

print("sample with replacement")
print(sample(sample30, replace=T))

readline()
print(sort(sample30)) 

readline()
print(sort(sample(sample30, replace=T)))

readline()

m <- 1000
res <- numeric(m)
for (i in 1:m) res[i] <- mean(sample(sample30, replace=T))
print("sd of the mean estimate")
print(sqrt(var(res)))                          # original size is just 30 obs.

readline()

# theoretical value s/sqrt(n)
print("theoretical result")
print(1/sqrt(30))

readline()

## distribution of means

print("distirbution of means of resamples")
hist(res, freq=FALSE)
lines(density(res))

readline()

# c.i.

print("c.i.")

print(quantile(res, c(0.025, 0.975)))

readline()


#### curve fitting

x <- seq(1,100,,200)
y <- c(x[x<70], x[x >= 70]^1.07) + rnorm(200,sd=7)

plot(x,y)
lines(lowess(x,y))
lines(x,fitted(lm(y ~ x)))

readline()

## bootstrap

indice <- 1:200
set.seed(33); m <- 3
indice.b <- vector("list",m)
lm.b <- vector("list",m)
loess.b <- vector("list",m)

for (i in 1:m) {
indice.b[[i]] <- sample(indice, replace=T)
lm.b[[i]] <- lm(y[indice.b[[i]]] ~ x[indice.b[[i]]])
loess.b[[i]] <- lowess(x[indice.b[[i]]],y[indice.b[[i]]])
}

## plot
#lm
plot(x,y)
for (i in 1:m) {
ordem <- order( x[ indice.b[[i]] ] )
lines( x[ indice.b[[i]] ][ordem], fitted( lm.b[[i]] )[ordem] )
}

readline()

plot(x,y, type="n")
#loess
for (i in 1:m) {
ordem <- order( x[ indice.b[[i]] ] )
#lines( x[ indice.b[[i]] ][ordem], predict( loess.b[[i]] )[ordem] )
lines( loess.b[[i]] )[ordem]
}