rm(list=ls())

pdf(file="teoriaassintotica.pdf",width=9,height=6)

set.seed(7594659)
as    = c(1,1,9)
bs    = c(9,1,1)
par(mfrow=c(1,3))
for (i in 1:3){
	mean  = as[i]/(as[i]+bs[i])
	var   = as[i]*bs[i]/((as[i]+bs[i])^2*(as[i]+bs[i]+1))
  error = (rbeta(100000,as[i],bs[i])-mean)/sqrt(var)
  plot(density(error),xlab="",main="")
  title(paste("Erro ",i,sep=""))
}

M     = 1000
ns    = c(10,50,100)
beta  = 2.0
coef  = array(0,c(3,3,M))
for (i in 1:3){
  a = as[i]
  b = bs[i]
  for (j in 1:3){
  	  n = ns[j]
  	  for (k in 1:M){
      x        = rnorm(n)
      error    = 6*rbeta(n,a,b)-3
      y        = beta*x + error
      coef[i,j,k] = lm(y~x-1)$coef
    }
  }
}


par(mfrow=c(3,3))
for (j in 1:3)
  for (i in 1:3){
    hist(coef[i,j,],xlab="",main="",prob=TRUE,
         xlim=range(coef))
    title(paste("n=",ns[j],"\nError ",i,sep=""))
    abline(v=beta,col=2)    
}



set.seed(987553)
error1 = rbeta(100000,1,9)
error2 = rbeta(100000,10,10)
error3 = rbeta(100000,9,1)
ind = sample(1:3,size=100000,replace=TRUE,prob=rep(1/3,3))
error = error1
error[ind==2]=error2[ind==2]
error[ind==3]=error3[ind==3]
error = 5*(error-0.5)
par(mfrow=c(2,2))
plot(density(error),xlab="",main="Error 4")

M     = 1000
ns    = c(10,50,100)
beta  = 2.0
coef  = matrix(0,M,3)
for (j in 1:3){
  n = ns[j]
  for (k in 1:M){
    x = rnorm(n)
    error1 = rbeta(n,1,9)
    error2 = rbeta(n,10,10)
    error3 = rbeta(n,9,1)
    ind = sample(1:3,size=n,replace=TRUE,prob=rep(1/3,3))
    error = error1
    error[ind==2]=error2[ind==2]
    error[ind==3]=error3[ind==3]
    error = 5*(error-0.5)
    y        = beta*x + error
    coef[k,j] = lm(y~x-1)$coef
  }
}

for (j in 1:3){
  hist(coef[,j],xlab="",main="",prob=TRUE,xlim=c(0,4))
  title(paste("n=",ns[j],"\nError 4",sep=""))
  abline(v=beta,col=2)    
}


dev.off()