Ejercicio 1 : Automatización de Análisis
El archivo DatosTarea2.zip contiene 332 archivos .csv con datos de monitoreo de contaminación del aire por partículas finas (PM) en 332 ubicaciones de cierto país. Cada archivo contiene datos un monitor y el número de identificación de cada monitor está incluido en el nombre del archivo. Por ejemplo, los datos del monitor 200 están en el archivo “200.csv”. Cada archivo contiene tres variables:
- Date: la fecha de la observación en formato AAAA-MM-DD (año-mes-día).
- sulfate: el nivel de PM de sulfato en el aire en esa fecha (medido en microgramos por metro cúbico).
- nitrate: el nivel de PM de nitrato en el aire en esa fecha (medido en microgramos por metro cúbico).
En cada archivo, notará que hay muchos días en los que faltan sulfato o nitrato (o ambos) (codificados como NA). Esto es común con los datos de monitoreo de la contaminación del aire. Todos los días nos piden informes para una lista específica de monitores y para una determinada variable. Se requiere construir proceso automatizado que lea el directorio donde se encuentran todos los datos de los monitores cargue los datos y entregue el análisis solicitado.
Parte 1: Creación de función
Primero describimos la presencia de NAs en la base de datos como se muestra a continuación en la siguiente tabla:
# librerias necesarias
library(tidyverse)
library(naniar)
library(formattable)
library(dplyr)
library(lubridate)
library(knitr)
library(graphics)
archivos <- list.files("DatosTarea2")
archivos <- paste0("DatosTarea2\\", archivos)
# Método del for
pm.for <- tibble()# Creamos un tibble vacío
for(i in seq_along(archivos)){
pm.for <- pm.for %>%
bind_rows(read.csv(file = archivos[i]))
}
knitr::kable(miss_var_summary(pm.for),
caption = 'Tabla resumen de los NA por Variables')| variable | n_miss | pct_miss |
|---|---|---|
| nitrate | 657738 | 85.18962 |
| sulfate | 653304 | 84.61533 |
| Date | 0 | 0.00000 |
| ID | 0 | 0.00000 |
Se observa que en base al contexto del ejercicio, el 85.19% equivales a valores del contaminante nitrato y el 84.62% correspondiente al sulfato, no presentan registros (NAs) en la base completa de los monitores.
### Función Generadora de Resumen -------------------------------------------
Resumen_Monitor <- function(x,y){
if(x >= 1 & x <= 332){
if(y =="sulfato" | y == "nitrato"){
if(y == "sulfato"){
df_resumen <- pm.for %>%
separate(Date, c("anio", "mes", "dia")) %>%
select(-dia) %>%
filter(ID == x) %>%
group_by(ID = x,anio,mes) %>%
summarise(Porcentaje_Casos_Completos = percent(n_complete(sulfate)/n(),1),
Promedio_Mensual = mean (sulfate, na.rm = TRUE)) %>%
as.data.frame()
df_resumen
}else{
df_resumen <- pm.for %>%
separate(Date, c("anio", "mes", "dia")) %>%
select(-dia) %>%
filter(ID == x) %>%
group_by(ID = x,anio,mes) %>%
summarise(Porcentaje_Casos_Completos = percent(n_complete(nitrate)/n(),1),
Promedio_Mensual = mean (nitrate, na.rm = TRUE)) %>%
as.data.frame()
df_resumen
}
}else{
print("ingrese correctamente si desea el reporte de sulfato o nitrato")
}
}else{
"ingrese correctamente el monitor a reportar (entre 1 a 332) "
}
}- Pruebe su función con el monitor 20 y con la variable nitrato. Presente en su informe el head() y tail() del marco de datos resultante.
NOTA: Se decicio dejar la variable ID del monitor para poder discriminar al aplicar la funcion con distintos argumentos en caso que corresponda.
knitr::kable(head(Resumen_Monitor(20,"nitrato")),
caption = 'Tabla: head() del Monitor 20, analizando nitrato')| ID | anio | mes | Porcentaje_Casos_Completos | Promedio_Mensual |
|---|---|---|---|---|
| 20 | 2002 | 01 | 0.0% | NaN |
| 20 | 2002 | 02 | 10.7% | 0.7336667 |
| 20 | 2002 | 03 | 12.9% | 0.9670000 |
| 20 | 2002 | 04 | 10.0% | 0.4460000 |
| 20 | 2002 | 05 | 12.9% | 0.6197500 |
| 20 | 2002 | 06 | 13.3% | 0.4122500 |
knitr::kable(tail(Resumen_Monitor(20,"nitrato")),
caption = 'Tabla: tail() del Monitor 20, analizando nitrato')| ID | anio | mes | Porcentaje_Casos_Completos | Promedio_Mensual | |
|---|---|---|---|---|---|
| 43 | 20 | 2005 | 07 | 0.0% | NaN |
| 44 | 20 | 2005 | 08 | 0.0% | NaN |
| 45 | 20 | 2005 | 09 | 0.0% | NaN |
| 46 | 20 | 2005 | 10 | 0.0% | NaN |
| 47 | 20 | 2005 | 11 | 0.0% | NaN |
| 48 | 20 | 2005 | 12 | 0.0% | NaN |
- Pruebe su función con el monitor 250 y con la variable sulfato. Presente en su informe el head() y tail() del marco de datos resultante.
knitr::kable(head(Resumen_Monitor(250,"sulfato")),
caption = 'Tabla: tail() del Monitor 250, analizando sulfato')| ID | anio | mes | Porcentaje_Casos_Completos | Promedio_Mensual |
|---|---|---|---|---|
| 250 | 2002 | 01 | 0.0% | NaN |
| 250 | 2002 | 02 | 0.0% | NaN |
| 250 | 2002 | 03 | 0.0% | NaN |
| 250 | 2002 | 04 | 13.3% | 3.8225 |
| 250 | 2002 | 05 | 16.1% | 4.6020 |
| 250 | 2002 | 06 | 13.3% | 9.0475 |
knitr::kable(tail(Resumen_Monitor(250,"sulfato")),
caption = 'Tabla: tail() del Monitor 250, analizando sulfato')| ID | anio | mes | Porcentaje_Casos_Completos | Promedio_Mensual | |
|---|---|---|---|---|---|
| 43 | 250 | 2005 | 07 | 12.9% | 6.8775 |
| 44 | 250 | 2005 | 08 | 16.1% | 6.8400 |
| 45 | 250 | 2005 | 09 | 0.0% | NaN |
| 46 | 250 | 2005 | 10 | 0.0% | NaN |
| 47 | 250 | 2005 | 11 | 0.0% | NaN |
| 48 | 250 | 2005 | 12 | 0.0% | NaN |
En ambos casos se puede observar la presencia de los NAs dentro de la misma, y por tanto para los próximos análisis no se consideraran.
Parte 2: Aplicación de la función y generalización.
Se realiza una función que en base a argumentos como: directorio, monitor y contaminante, retorne una base de datos que se encuenta en cierto directorio, con el promedio de los datos completos y promedio, análisis agregado por año y mes, que se llamará unir_monitores como a continuación se presenta:
unir_monitores <- function(directorio,monitor,contaminante){
archivos <- list.files(directorio)
archivos <- paste0(directorio,"\\", archivos)
pm.for <- tibble()
for(i in seq_along(archivos)){
pm.for <- pm.for %>%
bind_rows(read.csv(file = archivos[i]))
}
if(min(monitor) >= 1 &
max(monitor)<= 332 &
(contaminante == "nitrato" | contaminante == "sulfato")){
monitor <- as.integer(monitor)
i = 2
df_final <- Resumen_Monitor(monitor[1],contaminante)
while(i <= length(monitor)){
df_final <- df_final %>%
bind_rows(Resumen_Monitor(monitor[i],contaminante))
i = i+1
}
df_final
}else{
print("Favor de ingresar correctamente los argumento \n")
print("El directorio se refiere la carpeta donde se encuentran los archivos \n")
print("monitor se refiere al número de este, que va de 1 a 332 \n")
print("contaminante puede ser: nitrato o sulfato")
}
}Luego se procede a guardar en un objeto “Resumen”, con el fin de generar un gráfico de linea que represente la evolución de los promedios agregado por mes de los monitores y contaminantes escojidos para analizar.
Resumen1 <- unir_monitores("DatosTarea2",1:5,"nitrato")
Plot_resumen1 <- Resumen1 %>%
group_by(mes) %>%
summarise(Promedio_Mensual = mean (Promedio_Mensual, na.rm = TRUE)) %>%
as.data.frame() %>%
select(mes,Promedio_Mensual) %>%
plot(type = "l",
main = " Evolución del promedio mensual agregada del nitrato \n de los monitores 1 al 5",
ylab = "Promedio de nitrato en el aire",
xlab = "Meses",
lwd = 2,
col = "blue")
Plot_resumen1## NULLinicio <- Sys.time()
Resumen2 <- unir_monitores("DatosTarea2",1:332,"nitrato")
Plot_resumen2 <- Resumen2 %>%
group_by(mes) %>%
summarise(Promedio_Mensual = mean (Promedio_Mensual, na.rm = TRUE)) %>%
as.data.frame() %>%
select(mes,Promedio_Mensual) %>%
plot(type = "l",
main = " Evol. del prom. mensual agregada del nitrato \n de los monitores 1 al 332",
ylab = "Promedio de nitrato en el aire",
xlab = "Meses",
lwd = 2,
col = "Red")
Plot_resumen2
Sys.time() - inicio # Demoró 45.29757 minsinicio <- Sys.time()
Resumen2 <- unir_monitores("DatosTarea2",1:332,"nitrato")
Plot_resumen2 <- Resumen2 %>%
group_by(mes) %>%
summarise(Promedio_Mensual = mean (Promedio_Mensual, na.rm = TRUE)) %>%
as.data.frame() %>%
select(mes,Promedio_Mensual) %>%
plot(type = "l",
main = " Evol. del prom. mensual agregada del nitrato \n de los monitores 1 al 332",
ylab = "Promedio de nitrato en el aire",
xlab = "Meses",
lwd = 2,
col = "Red")
Plot_resumen2
Sys.time() - inicio # Demoró 45.29757 mins
Parte 3: Aplicación con funciones apply y map.
unir_monitores_map <- function(directorio,monitor,contaminante){
archivos <- list.files("DatosTarea2")
archivos <- paste0(directorio,"\\", archivos)
pm.map <- map_dfr(.x = archivos, .f = read.csv)# me da un dataframe
if(min(monitor) >= 1 &
max(monitor)<= 332 &
(contaminante == "nitrato" | contaminante == "sulfato")){
monitor <- as.integer(monitor)
i = 2
df_final <- Resumen_Monitor(monitor[1],contaminante)
while(i <= length(monitor)){
df_final <- df_final %>%
bind_rows(Resumen_Monitor(monitor[i],contaminante))
i = i+1
}
df_final
}else{
print("Favor de ingresar correctamente los argumento \n")
print("El directorio se refiere la carpeta donde se encuentran los archivos \n")
print("monitor se refiere al número de este, que va de 1 a 332 \n")
print("contaminante puede ser: nitrato o sulfato")
}
}inicio <- Sys.time()
Resumen_map <- unir_monitores_map("DatosTarea2",1:332,"nitrato")
Sys.time() - inicio # Demoró 32 mins.Clara mente se observa una reducción de tiempo de 49 min a 32 min usando funciones map.
Ejercicio 2 : Teorema Central del límite.
Parte 1: Simulación de una Función Exponencial
Genere 1.000.000 observaciones de la variable aleatoria exponencial con parámetro λ=0.1, guardela en un vector llamado población, puede usar la función rexp(n,rate) que genera n variables aleatorias con distribución exponencial de pámetro λ = rate. Construya un histograma de la distribución (si usa la función hist(), pruebe con nclass=50). Explique la forma de su histograma. ¿Es consistente con la forma de la distribución exponencial?. Calcule la media y la varianza, ¿son similares a la esperanza y la varianza de la distribución?. Si son consistentes, guárdelas en un objeto llamado mu y sigma2, respectivamente. Use como semilla el número 3308.
poblacion <- rexp(1000000,0.1)
hist(poblacion, nclass=50, probability = TRUE)
Al observar el histograma, la forma si es consistente con la familia de distribuciones exponenciales, ya que su memoria de cálculo se puede resumir de la siguiente manera, a medida que aumentan los valores de x, el valor de y va disminuyendo exponencialmente y concuerda con el gráfico.
Al comparar la media = 10.01016 y su varianza = 100.0847, en tal población generada, con su esperanza = 10 y varianza que es = 100, se puede decir que son aproximadamente iguales, lo que es lógico, ya que se simulón sin errores una población con tal parámetro.
Por lo tanto, como son consistentes \(\mu_2 = 10.01016\) y \(\sigma_2^2 = 100.0847\).
mu2 <- mean(poblacion)
mu2## [1] 10.005sigma2 <- var(poblacion)
sigma2## [1] 99.82732Parte 2: Comparación de iteraciones con muestra normal estandarizada.
Tome 1000 muestras de tamaño n=30 de su vector poblacion. Calcule el promedio de cada muestra. Calcule un vector con los 1000 valores Zn de cada muestra. Construya un histograma de los valores Zn. Compare con la distribución Normal(0,1) ajustando una curva normal a su gráfico o realizando una gráfico con qqnorm().
Z_n <- c()
for(i in 1:1000){
set.seed(3308+i)
muestra <- sample(poblacion,30, replace = T)
Z_n[i] <- (mean(muestra) - mu2)/sqrt(sigma2/30)
}
set.seed(3308)
hist(Z_n, nclass = 50, probability = TRUE,
main = "Curva densidad",
ylab = "Densidad")
lines(density(rnorm(30,0,1)), lwd = 2, col = 'red')
Claramente en el histograma se observa que las densidades de la normal estandarizada con la distribución del vector generado de una población exponencial se ajustan bien entre ellas, concluyendo así que las medias de las poblaciones generadas, ante muestras grandes, distribuyen normal con media 0 y varianza 1.
Parte 3: Creación de una Función(datos, tamaño de muestra, Iteraciones)
Cree una función que tome s muestras de tamaño n de su vector población, y devuelva un vector con los valores Zn de cada muestra y un histograma de los valores Zn con la distribución normal ajustada o un gráfico de qqnorm.
Distr_exp <- function(datos, n, s){
Z_n <- c()
mu3 <- mean(datos)
sigma3 <- var(datos)
for(i in 1:s){
set.seed(3308+i)
muestra <- sample(datos,n, replace = T)
Z_n[i] <- (mean(muestra) - mu3)/sqrt(sigma3/n)
}
mu4 <- mean(Z_n)
sigma4 <- var(Z_n)
hist(Z_n, nclass = 50, probability = TRUE,
main = "Curva densidad",
ylab = "Densidad")
lines(density(rnorm(n,mu4,sigma4)), lwd = 2, col = 'red')
}Parte 4: Puesta en práctica 1 de la función generada en la parte 3.
Ejecute su función para s = 1000 muestras de tamaño n = 5, 50, 100 y 1000 de su vector población. Compare los histogramas, ¿Qué diferencia se nota cuando aumenta el tamaño de la muestra en una población?
par(mfrow=c(2,2))
Distr_exp(poblacion,5,1000)
Distr_exp(poblacion,50,1000)
Distr_exp(poblacion,100,1000)
Distr_exp(poblacion,1000,1000)
Al ver que el tamaño de muestra va en aumento, la distribución de los promedios de la población exponencial se ajusta de mejor manera a una normal de media 0 y varianza 1.
Parte 5: Puesta en práctica 2 de la función generada en la parte 3.
Ejecute su función para s = 1000 muestras de tamaño n = 5, 50, 100 y 1000 en una población cuya distribución es exponencial con parámetro λ=10
set.seed(3308)
Exponencial <- rexp(1000000,10)
par(mfrow=c(2,2))
Distr_exp(Exponencial,5,1000)
Distr_exp(Exponencial,50,1000)
Distr_exp(Exponencial,100,1000)
Distr_exp(Exponencial,1000,1000)