VISUALIZANDO EL TEOREMA CENTRAL DEL LÍMITE CON R

 

Introducción

En este post, vamos a echarle un vistazo al TCL (teorema del límite central).

El TCL nos dice que las medias muestrales siguen, en el límite, una distribución normal. Pero - y esto es importante - no nos garantiza velocidad de convergencia.

Veremos dos casos:

• El caso “bueno”(t-student), donde la distribución de partida es simétrica y unimodal, y por lo tanto se espera una convergencia rápida.
• El caso malo (exponencial), una distribución muy sesgada, donde la convergencia de las medias muestrales es mucho más lenta.

Como siempre, está todo parametrizado y comentado de forma que - al menos esa es la intención - sea fácil modificar el código de acuerdo con tus necesidades.

Carga de librerías, parámetros globales del programa

Comenzaremos definiendo los parámetros globales. La tasa de la exponencial, grados de libertad de la t-student, número de muestras a tomar…

# Necesitamos esta librería para el test de normalidad (lilliefors test)
library(nortest)

# para que todo tengamos los mismos resultados
set.seed(123)

#número de muestras que tomaremos para cada tamaño muestral
n_samples_per_size<-1000

lambda<- 10 # parametro de la exponencial
n<-5 # grados de libertad de la t-student

Muestras de una t-student

Vamos a generar una buena cantidad de muestras, n_samples_per_size, para cada tamaño diferente. Calcularemos la media de cada una de ellas y la guardaremos en nuestro vector de medias muestrales. Queremos saber qué distribución tiene este vector.

En este caso, la distribución de partida es una t-student con pocos grados de libertad (hasta df=20-30 la t-student y la normal no tienen demasiado parecido). Esto es así para que se vea que hay convergencia rápida incluso en un caso no tan óptimo como podría ser df=30.

Esta es la t-student, para distintos grados de libertad, comparada con la normal. Puede visualizarse como una “normal con colas más pesadas y algo más plana en el centro”.

x <- seq(-3, 3, length=1000)
hx <- dnorm(x)

degf <- c(1, 2, 5)
colors <- c("red", "blue", "darkgreen", "black")
labels <- c("df=1", "df=2", "df=5", "normal")

plot(x, hx, type="l", lty=1, xlab="valor de x",
     ylab="densidad", main="comparando t-student",lwd=2)

for (i in 1:length(degf)){
    lines(x, dt(x,degf[i]), lwd=1, col=colors[i])
}

legend("topleft", inset=.05, title="Distributions",
       labels, lwd=1, lty=c(1, 1, 1, 1), col=colors)

Una vez completadas todas las muestras, veremos con un test de normalidad (Lilliefors) si realmente sus medias tienen una distribución normal y la representaremos graficamente junto con el p-valor asociado al test de normalidad.

for (sample_size in c(5,10,20,30))
{
avg<-rep(0,n_samples_per_size) # vector de las medias muestrales 
sdev<-rep(0,n_samples_per_size)


  for(i in 1:n_samples_per_size)
        {
        x<-rt(sample_size,n)
        avg[i]<-mean(x)
        sdev[i]<-sd(x)
  } 
  pvalor<-round(lillie.test(avg)$p.value,3)
  titulo<-paste("n=",sample_size, "media=", round(mean(avg),2), "sd=",
                round(sd(avg),2),"pvalor=",pvalor,"\n",sep=" ")
  hist(avg,main=titulo,col="red")
  }  

Como comentábamos antes, la distribución de la media muestral es muy parecida a la normal y ajusta muy bien desde un principio. Esto se debe a la forma “acampanada” de la t-student.

Veamos que pasa en el caso malo.

Muestras de una exponencial

Ahora vamos a repetetir el experimento con una distribución muy sesgada: la exponencial. Veremos que, aunque es cierto que converge a una normal, lo hace de una forma mucho más lenta. Tendremos que ir aumentando el tamaño muestral mucho más de lo esperado para conseguir una normalidad decente.

Esta es la curva de densidad de la distribución exponencial:

#exponencial original
titulo<-paste("exp(",lambda,")",sep="")
curve(dexp(x, rate = lambda), from=0, to=2, col='blue',main=titulo,ylab="")

Y, finalmente, el código para la exponencial, que es lo mismo que hacíamos arriba pero para esta nueva distribución.

Recordar que la media de la exp(10) es precisamente 0.1, que es el valor que nos sale en cada instancia.

sample_sizes_list <- c(5,10,30,50,75,100,250)
#par(mfrow=c(2,2))

for (sample_size in sample_sizes_list)
{
avg<-rep(0,n_samples_per_size) # vector de las medias muestrales 
  
  for(i in 1:n_samples_per_size)
        {
        x<-rexp(sample_size,lambda)
        avg[i]<-mean(x)
  } 
  pvalor<-round(lillie.test(avg)$p.value,4) #test normalidad para la media muestral
  titulo<-paste("n=",sample_size,"media=",round(mean(avg),2),"pvalor=", pvalor, sep=" ")
  hist(avg,main=titulo,col="blue")
}        

Se ve que le cuesta alcanzar la normalidad. Hasta el valor n=75 tiene demasiada asimetría para ser considerada como tal.

RESUMEN

En este post, hemos creado una herramienta en R para visualizar la convergencia a la normalidad de las medias muestrales - de acuerdo con el TCL (teorema central del límite).

Hemos visto como, para distrubuciones acampanadas (parecidas a la normal) tenemos una convergencia muy rápida, ya efectiva antes de n=30.

Sin embargo para distribuciones sesgadas la convergencia es mucho más lenta, teniéndonos que ir a valores cercanos a 100 para conseguir normalidad.