Teorema Central del Límite: Métodos Monte Carlo en R

In [1]:
# Muestra la fecha de la última ejecución del proyecto
print(paste0("Última ejecución: ", date()))
[1] "Última ejecución: Mon Nov 25 22:39:01 2019"

Contenido

  1. Introducción
  2. TCL
  3. Muestreo Aleatorio
    1. Ejercicio 1
  4. S. Monte Carlo
    1. Ejercicio 1
  5. Coméntanos










Introducción

Por muchos considerado el resultado más importante de las matemáticas, el teorema central del límite, en la teoría de la probabilidad, establece a la distribución Normal como la distribución a la que la media de casi cualquier conjunto de variables independientes y generadas aleatoriamente converge rápidamente. El teorema central del límite explica por qué la distribución Normal surge con tanta frecuencia y por qué generalmente es una excelente aproximación de la media de una recopilación de datos (a veces con tan sólo 10 variables).


Encyclopedia Britannica. (2019).

Central Limit Theorem.


En este proyecto desarrollaremos ejercicios para representar, ayudar a comprender y aplicar el teorema central del límite realizando experimentos Monte Carlo (aleatorios).
Todos los desarrollos de Pakin Ciencia de Datos son dinámicos por lo que periódicamente se añaden nuevos ejercicios a este proyecto. Al final de la página puedes comentarnos e incluso solicitar o sugerir el desarrollo de algún ejercicio o aplicación en particular.

Regresar al menú


Teorema Central del Límite

Sea $X_{1}, X_{2},...$ una sucesión infinita de variables aleatorias independientes e idénticamente distribuidas, con media $\mu$ y varianza finita $\sigma^2$. Entonces la función de distribución de la variable aleatoria:



$$ Z_{n} = \frac{(X_{1} +...+X_{n}) - n\mu}{\sqrt{n\sigma^2}}$$


tiende a la función de distribución normal estándar cuando n tiende a infinito.


Rincón, L. (2014), Introducción a la Probabilidad, México CDMX, Las prensas de ciencias.


Aquí puedes descargar el libro:

Introducción a la Probabilidad

En palabras más simples, el TCL nos dice que si tomamos muestras lo suficientemente grandes (>=30) de una cierta población, la media de esas muestras seguirá una distribución normal sin importar que distribución siga la población.

Regresar al menú


Muestreo Aleatorio

Ejercicio 1. Importación y Muestreo Aleatorio de Datos de Salud

Vamos a importar un conjunto de datos correspondientes a una encuesta de salud realizada en México por el Instituto Nacional de Salud Pública. Sólo importaremos el campo que contiene el peso (kg) de las personas encuestadas.

In [2]:
# Importamos los datos directamente del URL en donde se encuentran almacenados
datos <- read.csv("http://insp.mx/encuestoteca/Encuestas/ENSANUT2012_Nutricion/D%29Sobrepeso%20y%20obesidad_adulta_N2012%20%28CSV%29.csv",
                 header = TRUE)[ ,c('peso')]

# Eliminamos datos faltantes
datos <- na.omit((datos))

# Redondeamos a 2 decimales
datos <- round(datos, 2)

# exploramos el tamaño de la población
print("Cantidad de datos")

length(datos)

print("Observemos los datos")

head(datos,30)

print("Media de los pesos")

mp = round(mean(datos), 2)
mp
[1] "Cantidad de datos"
38267
[1] "Observemos los datos"
  1. 75
  2. 62.3
  3. 60.2
  4. 87.4
  5. 54.6
  6. 52.65
  7. 57.6
  8. 69.25
  9. 55.4
  10. 66.8
  11. 67.9
  12. 78.8
  13. 90.8
  14. 79.2
  15. 65.5
  16. 76.7
  17. 78.15
  18. 69
  19. 87.8
  20. 72.35
  21. 69.6
  22. 45.8
  23. 70.55
  24. 71.8
  25. 81.7
  26. 90.7
  27. 49.68
  28. 53.4
  29. 83.4
  30. 77.25
[1] "Media de los pesos"
70.86

Observemos la distribución:

In [3]:
hist(datos, main = "Distribución de pesos (encuesta)", xlab = "Peso (kg)", ylab = "Frecuencia", col = "green",
    border = "blue")

Hagamos experimentos Monte Carlo tomando n muestras aleatorias de tamaño m:

In [4]:
# Declaramos 2 variables n y m

# Número de muestras
n <- 10

# Tamaño de la muestra
m <- 10

# Creamos un vector de longitud m con entradas numéricas en donde almacenaremos las medias de los pesos
Medias <- vector("numeric", length = m)

# Creamos un ciclo para tomar m muestras de tamaño n y calcular su media
for(i in 1:m){

    Medias[i] <- mean(sample(datos, size = n, replace = FALSE))
}

print("Observemos las Medias")

head(round(Medias, 2),30)

print("Media de los pesos")

mp

print("Media de las muestras")

round(mean(Medias), 2)
[1] "Observemos las Medias"
  1. 66.61
  2. 73.54
  3. 65.68
  4. 69.73
  5. 71.19
  6. 71.05
  7. 65.48
  8. 75.2
  9. 69.66
  10. 70
[1] "Media de los pesos"
70.86
[1] "Media de las muestras"
69.81
In [5]:
# Ajustamos el tamaño de las gráficas
options(repr.plot.width=6, repr.plot.height=5)

Graficamos el histograma de las medias obtenidas:

In [6]:
hist(Medias, main = "Distribución de Medias", col= "red", ylab = "Frecuencia")

Incrementemos el número de muestras:

In [7]:
n <- 50
m <- 10

Medias <- vector("numeric", length = m)

for(i in 1:m){
    
    Medias[i] <- mean(sample(datos, size = n, replace = FALSE))
}

print("Observemos las Medias")

head(round(Medias,2),30)

print("Media de los pesos")

mp

print("Media de las muestras")

round(mean(Medias), 2)
[1] "Observemos las Medias"
  1. 71.56
  2. 70.56
  3. 71.95
  4. 70.99
  5. 71.1
  6. 71.61
  7. 73.44
  8. 69.46
  9. 67.8
  10. 71.49
[1] "Media de los pesos"
70.86
[1] "Media de las muestras"
71
In [8]:
hist(Medias, main = "Distribución de Medias", col= "red", ylab = "Frecuencia")

Incrementemos el número y el tamaño de las muestras:

In [9]:
n <- 100
m <- 100

Medias <- vector("numeric", length = m)

for(i in 1:m){
    
    Medias[i] <- mean(sample(datos, size = n, replace = FALSE))
}

print("Observemos las Medias")

head(round(Medias,2),30)

print("Media de los pesos")

mp

print("Media de las muestras")

round(mean(Medias), 2)
[1] "Observemos las Medias"
  1. 68.13
  2. 71.15
  3. 71.19
  4. 72.03
  5. 69.96
  6. 71.21
  7. 68.37
  8. 70.4
  9. 72.27
  10. 70.79
  11. 70.05
  12. 69.51
  13. 70.36
  14. 71.8
  15. 68.24
  16. 69.77
  17. 70.29
  18. 71.41
  19. 71.64
  20. 72.49
  21. 67.61
  22. 69.3
  23. 72.25
  24. 72.27
  25. 69.57
  26. 71.24
  27. 71.27
  28. 71.98
  29. 70.36
  30. 71.53
[1] "Media de los pesos"
70.86
[1] "Media de las muestras"
70.71
In [10]:
hist(Medias, main = "Distribución de Medias", col= "red", ylab = "Frecuencia")

Observamos como el histograma anterior comienza a normalizarse.

Incrementemos drásticamente el número y el tamaño de las muestras:

In [11]:
n <- 20000
m <- 2000

Medias <- vector("numeric", length = m)

for(i in 1:m){

    Medias[i] <- mean(sample(datos, size = n, replace = FALSE))
}

print("Observemos las Medias")

head(round(Medias,2),30)

print("Media de los pesos")

mp

print("Media de las muestras")

round(mean(Medias), 2)
[1] "Observemos las Medias"
  1. 70.91
  2. 70.91
  3. 70.83
  4. 70.94
  5. 70.75
  6. 71.01
  7. 70.9
  8. 70.75
  9. 70.89
  10. 70.93
  11. 70.64
  12. 70.73
  13. 70.91
  14. 70.99
  15. 71.03
  16. 70.71
  17. 70.85
  18. 70.88
  19. 70.84
  20. 70.88
  21. 70.88
  22. 70.92
  23. 70.78
  24. 70.85
  25. 70.9
  26. 70.83
  27. 70.86
  28. 70.78
  29. 70.91
  30. 70.95
[1] "Media de los pesos"
70.86
[1] "Media de las muestras"
70.86
In [12]:
hist(Medias, main = "Distribución de Medias", col= "red", ylab = "Frecuencia")

  1. ¿Qué podemos decir acerca de esta última distribución?
  2. ¿Qué podemos decir de su media?
  3. ¿Qué podemos decir de su desviación estándar?

Regresar al menú


Sumulaciones Monte Carlo

Los métodos Monte Carlo son procedimientos no deterministas que pueden ser utilizados para aproximar y/o simular variables aleatorias, su nombre hace referencia al Casino de Monte Carlo en Mónaco por ser un ícono de los juegos de azar.

Ejercicio 1. Simulación Monte Carlo Lanzamientos de una Moneda

Con la creación de un sencillo ciclo for simularemos lanzamientos de una moneda (justa) en donde se representará la cara de la moneda con el número 1 y la cruz con el numero 2. Haremos m repeticiones de n lanzamientos cada repetición y nos fijaremos en el número de veces que aparece la cara (1) y promediaremos (media) este número de veces que a su vez acumularemos en un vector llamado Resultados.

Comencemos con 10 repeticiones de 10 lanzamientos cada repetición:

In [13]:
n <- 10
m <- 10

Resultados <- vector("numeric", length = m)

for(i in 1:m){
    
Resultados[i] <- mean(sample(1:2, size = n, replace = TRUE) == 1)

}

Resultados
  1. 0.5
  2. 0.4
  3. 0.4
  4. 0.7
  5. 0.3
  6. 0.6
  7. 0.6
  8. 0.3
  9. 0.4
  10. 0.4
In [14]:
hist(Resultados, main = "Distribución de Medias (Simulación Monte Carlo)",
     col= "chocolate2", ylab = "Frecuencia")

Incrementemos el número de lanzamiento por cada repetición:

In [15]:
n <- 100
m <- 10

Resultados <- vector("numeric", length = m)

for(i in 1:m){
    
Resultados[i] <- mean(sample(1:2, size = n, replace = TRUE) == 1)

}

Resultados
  1. 0.48
  2. 0.53
  3. 0.51
  4. 0.46
  5. 0.43
  6. 0.53
  7. 0.41
  8. 0.51
  9. 0.5
  10. 0.47
In [16]:
hist(Resultados, main = "Distribución de Medias (Simulación Monte Carlo)",
     col= "chocolate2", ylab = "Frecuencia")

Incrementemos el número de repeticiones:

In [17]:
n <- 100
m <- 100

Resultados <- vector("numeric", length = m)

for(i in 1:m){
    
Resultados[i] <- mean(sample(1:2, size = n, replace = TRUE) == 1)

}

head(Resultados, 20)
  1. 0.52
  2. 0.4
  3. 0.49
  4. 0.48
  5. 0.42
  6. 0.53
  7. 0.45
  8. 0.49
  9. 0.62
  10. 0.42
  11. 0.54
  12. 0.54
  13. 0.46
  14. 0.56
  15. 0.49
  16. 0.59
  17. 0.47
  18. 0.53
  19. 0.58
  20. 0.41
In [18]:
hist(Resultados, main = "Distribución de Medias (Simulación Monte Carlo)",
     col= "chocolate2", ylab = "Frecuencia")

Vemos ya como el histograma anterior comienza a mostrar una distribución conocida.

Incrementemos drásticamente el número de lanzamientos y el número repeticiones:

In [19]:
n <- 10000
m <- 10000

Resultados <- vector("numeric", length = m)

for(i in 1:m){
    
Resultados[i] <- mean(sample(1:2, size = n, replace = TRUE) == 1)

}

head(Resultados, 20)
  1. 0.5024
  2. 0.4992
  3. 0.5015
  4. 0.5057
  5. 0.5033
  6. 0.5018
  7. 0.5032
  8. 0.4998
  9. 0.4962
  10. 0.5015
  11. 0.513
  12. 0.4997
  13. 0.5077
  14. 0.4964
  15. 0.4965
  16. 0.4952
  17. 0.4985
  18. 0.516
  19. 0.4936
  20. 0.5008
In [20]:
hist(Resultados, main = "Distribución de Medias (Simulación Monte Carlo)",
     col= "chocolate2", ylab = "Frecuencia")

  1. ¿Qué podemos decir acerca de esta última distribución?
  2. ¿Qué podemos decir de su media?
  3. ¿Qué podemos decir de su desviación estándar?

Regresar al menú


Síguenos en Facebook


Pakin en Facebook

Esperamos con gusto tus comentarios y sugerencias para incrementar y mejorar este proyecto: