# Observando o comportamento amostral
# dos estimadores de minimos quadrados
# ordinarios
set.seed(23235)
n = 100
B = 1000
bt= c(2,0.5)
bhat = matrix(0,B,2)
for (i in 1:B){
  x = rnorm(n)
  y = bt[1]+bt[2]*x + rnorm(n,0,0.25)
  bhat[i,] = lm(y~x)$coef
}
par(mfrow=c(2,2))
hist(bhat[,1])
abline(v=bt[1],col=2,lwd=2)
hist(bhat[,2])
abline(v=bt[2],col=2,lwd=2)

sqrt(var(bhat[,1]))
sqrt(var(bhat[,2]))


# Na vida real, infelizmente, em geral so' 
# temos um conjunto de dados
n = 100
bt= c(2,0.5)
x = rnorm(n)
y = bt[1]+bt[2]*x + rnorm(n,0,0.25)
bhat = lm(y~x)$coef


M = 1000
bhat1 = matrix(0,M,2)
for (i in 1:M){
  ind = sample(1:n,size=n,replace=TRUE)
  bhat1[i,] = lm(y[ind]~x[ind])$coef
}


hist(bhat1[,1])
abline(v=bt[1],col=2,lwd=2)
abline(v=bhat[1],col=3,lwd=2)

hist(bhat1[,2])
abline(v=bt[2],col=2,lwd=2)
abline(v=bhat[2],col=3,lwd=2)

# Comparando os desvios-padroes
# dos estimadores teoricos e bootstrap
sqrt(var(bhat[,1]))
sqrt(var(bhat1[,1]))

sqrt(var(bhat[,2]))
sqrt(var(bhat1[,2]))


# Segundo exemplo: correlacao
n = 100
bt= c(2,0.75)
x = rnorm(n)
y = bt[1]+bt[2]*x + rnorm(n,0,0.25)

r = cor(x,y)
M = 1000
r1 = rep(0,M)
for (i in 1:M){
  ind = sample(1:n,size=n,replace=TRUE)
  r1[i] = cor(x[ind],y[ind])
}

par(mfrow=c(1,1))
hist(r1)
abline(v=r,col=2,lwd=2)