###############################################################################
#
#             Instrucciones para resolver el ejercicio de LGN y TCL en R 
#                           Probabilidades y Estadistica (M)
#
###############################################################################


Ejercicio 1: 

#a) primero generamos mil v.a. uniformes y luego hacemos el histograma correspondiente. 

xbarra<-runif(1000)
hist(xbarra)

#b) primero generamos dos v.a. uniformes o un vector uniforme bivariado, con ambas coordenadas independientes

aa<-runif(2)
xbarra2<-mean(aa)

#ahora repetimos 1000 veces
#con un for:
xbarra2<-rep(0,1000)
for (i in 1:1000)
{
aa<-runif(2)
xbarra2[i]<-mean(aa)
}

hist(xbarra2)

#o vectorialmente;
a1<-runif(2000)
#convertimos al vector en matriz
matriz2<-matrix(a1,nrow=2, ncol=1000)

#aplicamos el promedio por columnas
xbarra2<-apply(matriz2,2,mean)

hist(xbarra2)

#observar que si uno repite todo el ítem b) con otros dos mil datos, obtiene un histograma ligeramente diferente

#c) idem 5
matriz5<-matrix(runif(5000),nrow=5, ncol=1000)

#chequeamos
dim(matriz5)

xbarra5<-apply(matriz5,2,mean)

# Para comparar los 3 histogramas en un mismo gráfico, dividimos la pantalla en 6

par(mfrow=c(2,3))
hist(xbarra,probability=T)
hist(xbarra2,probability=T)
hist(xbarra5,probability=T)

#d) ídem 30
matriz30<-matrix(runif(30000),nrow=30, ncol=1000)
xbarra30<-apply(matriz30,2,mean)
hist(xbarra30,probability=T)


#e) ídem 500

matriz500<-matrix(runif(500000),nrow=500, ncol=1000)
xbarra500<-apply(matriz500,2,mean)
hist(xbarra500,probability=T)

# f) idem 1200
matriz1200<-matrix(runif(1200*1000),nrow=1000, ncol=1200)
xbarra1200<-apply(matriz1200,2,mean)
hist(xbarra1200,probability=T)

hist(xbarra,probability=T)
hist(xbarra2,probability=T)
hist(xbarra5,probability=T)
hist(xbarra30,probability=T)
hist(xbarra500,probability=T)
hist(xbarra1200,probability=T)

#Mejor aún, pongamos la misma escala en x para comparar

hist(xbarra,probability=T,xlim=c(0,1))
hist(xbarra2,probability=T,xlim=c(0,1))
hist(xbarra5,probability=T,xlim=c(0,1))
hist(xbarra30,probability=T,xlim=c(0,1))
hist(xbarra500,probability=T,xlim=c(0,1))
hist(xbarra1200,probability=T,xlim=c(0,1))

#notemos que queda distorsionado si no obligamos a que la escala en el eje y sea la misma
#hallemos el máximo valor para el eje y en el último grafico
par(mfrow=c(2,3))

MM<-max(hist(xbarra1200,plot=F)$intensities)

hist(xbarra,probability=T,xlim=c(0,1),ylim=c(0,MM))
hist(xbarra2,probability=T,xlim=c(0,1),ylim=c(0,MM))
hist(xbarra5,probability=T,xlim=c(0,1) ,ylim=c(0,MM))
hist(xbarra30,probability=T,xlim=c(0,1) ,ylim=c(0,MM))
hist(xbarra500,probability=T,xlim=c(0,1) ,ylim=c(0,MM))
hist(xbarra1200,probability=T,xlim=c(0,1) ,ylim=c(0,MM))

# o bien, hagamos los gráficos de estimaciones de la densidad los 6 grupos de datos en el mismo gráfico
par(mfrow=c(1,1))
plot(density(xbarra1200),xlim=c(0,1),lwd=2,main="Densidades estimadas para los promedios de n observaciones
N(0,1)")
lines(density(xbarra500),col=2,lwd=2)
lines(density(xbarra30),col=6,lwd=2)
lines(density(xbarra5),col=4,lwd=2)
lines(density(xbarra2),col=5,lwd=2)
lines(density(xbarra),col=3,lwd=2)
legend( 0.6, 30, c("n = 1", "n = 2","n = 5", "n = 30","n = 500","n = 1200"), col = c(3,5,4,6,2,1),lty = c(1,1, 1, 1,1,1), lwd =rep(2,6))

boxplot(xbarra,xbarra2,xbarra5,xbarra30,xbarra500,xbarra1200)

#medias muestrales
mean(xbarra)
mean(xbarra2)
mean(xbarra5)
mean(xbarra30)
mean(xbarra500)
mean(xbarra1200)

#varianzas muestrales

var(xbarra)
var(xbarra2)
var(xbarra5)
var(xbarra30)
var(xbarra500)
var(xbarra1200)

#g) Estandaricemos a una escala adecuada para poder apreciar el comportamiento aleatorio de xbarra

e1<-(xbarra-0.5)*sqrt(1)/sqrt(1/12)
e2<-(xbarra2-0.5)*sqrt(2)/sqrt(1/12)
e5<-(xbarra5-0.5)*sqrt(5)/sqrt(1/12)
e30<-(xbarra30-0.5)*sqrt(30)/sqrt(1/12)
e500<-(xbarra500-0.5)*sqrt(500)/sqrt(1/12)
e1200<-(xbarra1200-0.5)*sqrt(1200)/sqrt(1/12)

plot(density(e1200),lwd=2,main="Densidades estimadas para los promedios estandarizados segun TCL
de n observaciones N(0,1)",ylim=c(0,0.41))
lines(density(e500),col=2,lwd=2)
lines(density(e30),col=6,lwd=2)
lines(density(e5),col=4,lwd=2)
lines(density(e2),col=5,lwd=2)
lines(density(e1),col=3,lwd=2)
lines(seq(-4,4,by=0.01),dnorm(seq(-4,4,by=0.01)),col=7,lwd=2)
legend( 2.5, 0.3, c("n = 1", "n = 2","n = 5", "n = 30","n = 500","n = 1200","densidad N(0,1)"), col = c(3,5,4,6,2,1,7),lty = c(1,1,1,1,1,1,1), lwd =rep(2,7))



par(mfrow=c(2,3))

hist(e1,probability=T,xlim=c(-4,4),ylim=c(0,0.4))
lines(seq(-4,4,by=0.01),dnorm(seq(-4,4,by=0.01)),col=2,lwd=2)
hist(e2,probability=T,xlim=c(-4,4),ylim=c(0,0.4))
lines(seq(-4,4,by=0.01),dnorm(seq(-4,4,by=0.01)),col=2,lwd=2)
hist(e5,probability=T,xlim=c(-4,4) ,ylim=c(0,0.4))
lines(seq(-4,4,by=0.01),dnorm(seq(-4,4,by=0.01)),col=2,lwd=2)
hist(e30,probability=T,xlim=c(-4,4) ,ylim=c(0,0.4))
lines(seq(-4,4,by=0.01),dnorm(seq(-4,4,by=0.01)),col=2,lwd=2)
hist(e500,probability=T,xlim=c(-4,4) ,ylim=c(0,0.4))
lines(seq(-4,4,by=0.01),dnorm(seq(-4,4,by=0.01)),col=2,lwd=2)
hist(e1200,probability=T,xlim=c(-4,4) ,ylim=c(0,0.4))
lines(seq(-4,4,by=0.01),dnorm(seq(-4,4,by=0.01)),col=2,lwd=2)

boxplot(e1,e2,e5,e30,e500,e1200)


# h) Para generar observaciones con distribución Cauchy (no tiene esperanza ni varianza finita, es simétrica alrededor del cero, su densidad es f(x)= 1 / (pi*(1 + (x^2))  ), la instrucción es
#rcauchy(n)

matriz1<-matrix(rcauchy(1000),nrow=1, ncol=1000)
xbarra<-apply(matriz1,2,mean)
matriz2<-matrix(rcauchy(2000),nrow=2, ncol=1000)
xbarra2<-apply(matriz2,2,mean)
matriz5<-matrix(rcauchy(5000),nrow=5, ncol=1000)
xbarra5<-apply(matriz5,2,mean)
matriz30<-matrix(rcauchy(30000),nrow=30, ncol=1000)
xbarra30<-apply(matriz30,2,mean)
matriz500<-matrix(rcauchy(500000),nrow=500, ncol=1000)
xbarra500<-apply(matriz500,2,mean)
matriz1200<-matrix(rcauchy(1000*1200),nrow=1200, ncol=1000)
xbarra1200<-apply(matriz1200,2,mean)


#no se ve mucho en este gráfico, hay mucho peso en las colas.
mm<-max(c(density(xbarra)$y,density(xbarra2)$y,density(xbarra5)$y,density(xbarra30)$y,density(xbarra500)$y,density(xbarra1200)$y))

par(mfrow=c(1,1))
plot(density(xbarra1200),lwd=2,ylim=c(0,mm),main="Densidades estimadas para los promedios de n observaciones
Cauchy(0,1)")
lines(density(xbarra500),col=2,lwd=2)
lines(density(xbarra30),col=6,lwd=2)
lines(density(xbarra5),col=4,lwd=2)
lines(density(xbarra2),col=5,lwd=2)
lines(density(xbarra),col=7,lwd=2)
lines(seq(-4,4,by=0.01),dcauchy(seq(-4,4,by=0.01)),col=7,lwd=2)
legend( 100, 0.3, c("n = 1", "n = 2","n = 5", "n = 30","n = 500","n = 1200","densidad Cauchy(0,1)"), col = c(7,5,4,6,2,1,3),lty = c(1,1,1, 1, 1,1,1), lwd =rep(2,7))

#miremos lo mismo entre -4 y 4, solo para los n grandes

plot(density(xbarra1200),lwd=2,ylim=c(0,mm),xlim=c(-4,4),main="Densidades estimadas para los promedios de n observaciones
Cauchy(0,1)")
lines(density(xbarra500),col=2,lwd=2)
lines(density(xbarra30),col=6,lwd=2)
lines(seq(-4,4,by=0.01),dcauchy(seq(-4,4,by=0.01)),col=3,lwd=2)
legend(-4, 0.35, c( "n = 30","n = 500","n = 1200","densidad Cauchy(0,1)"), col = c(6,2,1,3),lty = c( 1, 1,1,1), lwd =rep(2,4),cex=0.8)

#medias muestrales
mean(xbarra)
mean(xbarra2)
mean(xbarra5)
mean(xbarra30)
mean(xbarra500)
mean(xbarra1200)

#varianzas muestrales

var(xbarra)
var(xbarra2)
var(xbarra5)
var(xbarra30)
var(xbarra500)
var(xbarra1200)

Ejercicio 2:

library(animation)
ani.options(interval = 0.03, nmax = 213)
quincunx()