Taller de Simulación en R - Unidad 2
Métodos y simulación estadística
Problema 1
Estimación del valor de \(\pi\)
La siguiente figura sugiere como estimar el valor de \(\pi\) con una simulación. En la figura, un circuito con un área igual a \(\pi/4\), está inscrito en un cuadrado cuya área es igual a 1. Se elige de forma aleatoria \(n\) puntos dentro del cuadrado . La probabilidad de que un punto esté dentro del círculo es igual a la fracción del área del cuadrado que abarca a este, la cual es \(\pi/4\). Por tanto, se puede estimar el valor de \(\pi/4\) al contar el número de puntos dentro del círculo, para obtener la estimación de \(\pi/4\) . De este último resultado se encontrar una aproximación para el valor de \(\pi\) .
Pasos sugeridos:
-
Genere \(n\) coordenadas \(x\): \(X_{1}\), . . . , \(X_{n}\). Utilice la distribución uniforme con valor mínimo de \(0\) y valor máximo de \(1\). La distribución uniforme genera variables aleatorias que tienen la misma probabilidad de venir de cualquier parte del intervalo \((0, 1)\).
-
Genere \(1000\) coordenadas \(y\) : \(Y_{1}, . . . , Y_{n}\), utilizando nuevamente la distribución uniforme con valor mínimo de \(0\) y valor máximo de \(1\).
-
Cada punto \((X_{i},Y_{i})\) se encuentra dentro del círculo si su distancia desde el centro \((0.5, 0.5)\) es menor a \(0.5\). Para cada par \((X_{i},Y_{i})\) determine si la distancia desde el centro es menor a \(0.5\). Esto último se puede realizar al calcular el valor \((X_{i}-0.5)^{2}+(Y_{i}-0.5)^{2}\), que es el cuadrado de la distancia, y al determinar si es menor que \(0.25\).
-
¿Cuántos de los puntos están dentro del círculo?
## [1] 801- ¿Cuál es su estimación de \(\pi\)?
## Se establece el número de puntos a analizar con una generación de muestras n=1000
n <- 1000
x <- runif(n, min = 0, max = 1)#Se crea una funciónn para calcular la distancia del par de coordenadas al centro del Circulo
dist= (x-0.5)^2+(y-0.5)^2
#Se crea una variable que cuente la distancia de las parejas que cumplen con la condicición que sean menor a 0.25
c= as.numeric(dist<=0.25)
La relación entre el área del círculo y el cuadrado es de π/4.El área del círculo lo componen los puntos que se encuentran dentro del mismo (pi/4), en cuanto al área del cuadrado lo componen todos los n puntos.
#Área del cuadrado = 1
#Área del circulo = pi*R^2 = pi/4
# La relación entre el círculo y el cuadrado es = (pi/4)/1
# conteo/n = pi/4, para lo cual se despeja pi y se determina que es aproximadamente:
conteo*4/n## [1] 3.204Simular con cantidades superiores de puntos
Teniendo en cuenta lo planteado anteriormente, se crea una función que contenga los pasos anteriores y genere el número de puntos que están dentro del área del cículo y la estimación de π, esto con el fin de demostrar la reducción del error y la similitud con el valor verdadero de pi.
distancia <- function(xi, yi) {
sqrt((xi - 0.5)^2 + (yi - 0.5)^2)
}
estPi <- function(puntos) {
x <- runif(n = puntos, min = 0, max = 1)
y <- runif(n = puntos, min = 0, max = 1)
dist_punt <- distancia(x, y)
evaluar <- as.numeric(dist_punt <= 0.5)
puntosCirculo <- sum(evaluar)
pi_estimado <- formatC(puntosCirculo * 4 / puntos, format = "f", digits = 3)
result <- list(Total_Puntos = puntos, Puntos_circulo = puntosCirculo, Estimacion_pi = pi_estimado)
return(result)
}
puntos <- c(1000, 10000, 100000)
resultado <- sapply(puntos, estPi)
resultado## [,1] [,2] [,3]
## Total_Puntos 1000 10000 1e+05
## Puntos_circulo 788 7852 78358
## Estimacion_pi "3.152" "3.141" "3.134"Bajo los escenarios de simulación anteriores se puede evidenciar cómo a medida que se aumentan la cantidad de puntos a evaluar, así mismo el valor estimado de pi es cada vez más cercano a su valor real.
Problema 2
Propiedades de los estimadores
La simulación ayuda a entender y validad las propiedades de los estimadores estadísticos como son. insesgadez, eficiencia y la consistencia principalmente. El siguiente problema permite evidenciar las principales características de un grupo de estimadores propuestos para la estimación de un parámetro asociado a un modelo de probabilidad.
Sean \(X_{1}\), \(X_{2}\), \(X_{3}\) y \(X_{4}\), una muestra aleatoria de tamaño \(n = 4\) cuya población la conforma una distribución exponencial con parámetro \(\theta\) desconocido. Determine las características de cada uno de los siguientes estimadores propuestos:
- \(\widehat{\theta_1} = \dfrac{X_{1} + X_{2}}{6} + \dfrac{X_3 + X_4}{3}\)
- \(\widehat{\theta_2} = \dfrac{(X_1 + 2 X_2 + 3 X_3 + 4 X_4)}{5}\)
- \(\widehat{\theta_3} = \dfrac{X_1 + X_2 + X_3 + X_4}{4}\)
- \(\widehat{\theta_4} = \dfrac{\min\{X_1,X_2,X_3,X_4\} + \max\{X_1,X_2,X_3,X_4\}}{2}\)
# Se establece una semilla para reproducibilidad
set.seed(12)
muestra<- 20 #muestra inicial
lambda <- .5 #Parámetro θ a encontrar, se supone un valor inicial.
# Creación de variables
x1 <- rexp(muestra, lambda)
x2 <- rexp(muestra, lambda)
x3 <- rexp(muestra, lambda)
x4 <- rexp(muestra, lambda)
# Se crea un dataframe consolidando las variables y el muestreo.
base <- data.frame(x1,x2,x3,x4)
#Crear las funciones para cada estimador:
est1 <- function(x1,x2,x3,x4){
est = (x1+x2)/6 + (x3+x4)/3
return(est)
}
est2 <- function(x1,x2,x3,x4){
est = (x1+2*x2+3*x3+4*x4)/5
return(est)
}
est3 <- function(x1,x2,x3,x4){
est = (x1+x2+x3+x4)/4
return(est)
}
est4 <- function(x1,x2,x3,x4){
est = (min(x1,x2,x3,x4)+ max(x1,x2,x3,x4))/2
return(est)
}
# Se evalua cada estimador con la base de muestreo creada
E1 <- apply(base, 1, function(row){
est1(row["x1"], row["x2"], row["x3"], row["x4"])
})
E2 <- apply(base, 1, function(row){
est2(row["x1"], row["x2"], row["x3"], row["x4"])
})
E3 <- apply(base, 1, function(row){
est3(row["x1"], row["x2"], row["x3"], row["x4"])
})
E4<- apply(base, 1, function(row){
est4(row["x1"], row["x2"], row["x3"], row["x4"])
})
# Se crea una función para simular los diferentes escenarios con las muestras, está función tendrá como parámetros el tamaño de la muestra y el valor lambda.
escenarios_est <- function(muestras, lambda){
#Creación de Muestras
x1 <- rexp(muestras, lambda)
x2 <- rexp(muestras, lambda)
x3 <- rexp(muestras, lambda)
x4 <- rexp(muestras, lambda)
base <- data.frame(x1,x2,x3,x4)
# Aplicar a la base cada una de las funciones creadas para los estimadores.
E1 <- apply(base, 1, function(row){
est1(row["x1"], row["x2"], row["x3"], row["x4"])
})
E2 <- apply(base, 1, function(row){
est2(row["x1"], row["x2"], row["x3"], row["x4"])
})
E3 <- apply(base, 1, function(row){
est3(row["x1"], row["x2"], row["x3"], row["x4"])
})
E4<- apply(base, 1, function(row){
est4(row["x1"], row["x2"], row["x3"], row["x4"])
})
datest <- data.frame(E1,E2,E3,E4)
return(datest)
}Generación de muestras
Muestra n = 20
M20 <-escenarios_est(20,.5)
mis_colores <- c("#061A40", "#0353A4", "#006DAA", "#003559")
boxplot(M20, col = mis_colores, outline = FALSE, xlab = "Estimadores", ylab = "Valores", main = "Boxplot de los Estimadores")
abline(h = 2, col = "red")
## Descriptive Statistics
## M20
## N: 20
##
## E1 E2 E3 E4
## ----------------- -------- -------- -------- --------
## Mean 2.03 4.25 1.97 2.33
## Std.Dev 1.18 2.67 0.96 1.26
## Min 0.28 0.61 0.26 0.32
## Q1 1.28 2.61 1.30 1.46
## Median 1.86 3.71 2.01 2.10
## Q3 2.55 5.11 2.31 2.86
## Max 5.31 12.01 4.12 6.16
## MAD 0.92 1.81 0.86 1.05
## IQR 1.23 2.33 0.97 1.34
## CV 0.58 0.63 0.49 0.54
## Skewness 1.18 1.44 0.63 1.17
## SE.Skewness 0.51 0.51 0.51 0.51
## Kurtosis 1.15 1.77 0.08 1.80
## N.Valid 20.00 20.00 20.00 20.00
## Pct.Valid 100.00 100.00 100.00 100.00Muestra n = 50
M50 <-escenarios_est(50,.5)
mis_colores <- c("#061A40", "#0353A4", "#006DAA", "#003559")
boxplot(M50, col = mis_colores, outline = FALSE, xlab = "Estimadores", ylab = "Valores", main = "Boxplot de los Estimadores")
abline(h = 2, col = "red")
## Descriptive Statistics
## M50
## N: 50
##
## E1 E2 E3 E4
## ----------------- -------- -------- -------- --------
## Mean 1.94 3.81 1.95 2.20
## Std.Dev 1.09 2.15 1.01 1.18
## Min 0.33 0.78 0.27 0.44
## Q1 1.32 2.54 1.33 1.34
## Median 1.61 3.29 1.68 1.89
## Q3 2.57 4.55 2.48 2.90
## Max 6.23 13.12 5.32 6.04
## MAD 0.87 1.48 0.94 0.92
## IQR 1.19 1.96 1.13 1.55
## CV 0.56 0.56 0.52 0.54
## Skewness 1.59 1.90 1.09 1.04
## SE.Skewness 0.34 0.34 0.34 0.34
## Kurtosis 3.46 5.31 1.11 0.77
## N.Valid 50.00 50.00 50.00 50.00
## Pct.Valid 100.00 100.00 100.00 100.00Muestra n = 100
M100 <-escenarios_est(100,.5)
mis_colores <- c("#061A40", "#0353A4", "#006DAA", "#003559")
boxplot(M100, col = mis_colores, outline = FALSE, xlab = "Estimadores", ylab = "Valores", main = "Boxplot de los Estimadores")
abline(h = 2, col = "red")
## Descriptive Statistics
## M100
## N: 100
##
## E1 E2 E3 E4
## ----------------- -------- -------- -------- --------
## Mean 2.02 4.07 2.06 2.43
## Std.Dev 1.09 2.24 1.07 1.39
## Min 0.43 0.76 0.36 0.50
## Q1 1.19 2.48 1.27 1.46
## Median 1.80 3.58 1.85 1.99
## Q3 2.62 5.48 2.74 3.18
## Max 4.85 10.22 4.86 7.48
## MAD 0.97 2.15 1.05 1.15
## IQR 1.41 2.93 1.46 1.68
## CV 0.54 0.55 0.52 0.57
## Skewness 0.75 0.80 0.65 1.05
## SE.Skewness 0.24 0.24 0.24 0.24
## Kurtosis -0.21 -0.07 -0.32 0.95
## N.Valid 100.00 100.00 100.00 100.00
## Pct.Valid 100.00 100.00 100.00 100.00Muestra n = 1000
M1000 <-escenarios_est(1000,.5)
mis_colores <- c("#061A40", "#0353A4", "#006DAA", "#003559")
boxplot(M1000, col = mis_colores, outline = FALSE, xlab = "Estimadores", ylab = "Valores", main = "Boxplot de los Estimadores")
abline(h = 2, col = "red")
## Descriptive Statistics
## M1000
## N: 1000
##
## E1 E2 E3 E4
## ----------------- --------- --------- --------- ---------
## Mean 2.01 4.02 2.01 2.36
## Std.Dev 1.07 2.24 1.04 1.31
## Min 0.20 0.42 0.18 0.21
## Q1 1.23 2.43 1.24 1.42
## Median 1.82 3.58 1.84 2.12
## Q3 2.60 5.14 2.61 2.98
## Max 6.60 14.22 6.22 10.27
## MAD 1.01 1.95 0.97 1.13
## IQR 1.37 2.70 1.37 1.55
## CV 0.53 0.56 0.51 0.55
## Skewness 1.07 1.21 0.99 1.26
## SE.Skewness 0.08 0.08 0.08 0.08
## Kurtosis 1.55 2.09 1.20 2.58
## N.Valid 1000.00 1000.00 1000.00 1000.00
## Pct.Valid 100.00 100.00 100.00 100.00Conclusiones
Después de analizar cada uno de los escenarios con diferentes valores de muestras se puede afirmar que los mejores estimadores son \(\widehat{\theta_1}\) y el \(\widehat{\theta_3}\) según las propiedades y lo que se puede evidenciar en su gráfico de caja:
El estimador 1 es el que más se acerca al valor esperado y esto se demostró a medida que se incrementaron las muestras, es decir es el que cuenta con mayor afinidad con la propiedad de Insesgadez.
El estimador 2 presenta una varianza alta y además de esto es bastante sesgado ya que el valor de la media estaba por encima del valor esperado (2), En cuanto a consistencia si se evidenció que a medida que incrementó la muestra, empezó acercarse al valor del parámetro.
El estimador 3 es el que presenta menor varianza, el cual se analizó a partir de la desviación estandar, es decir es el más eficiente. Ademas< tambien presenta la propiedad de insesgadez ya que la Media de la distribución del estimador se acercaba bastante al parámetro y esto se denotaba aún más a medida que aumentaga el tamaño de la muestra.
El estimador 4 no se encuentra tan lejos del valor esperado presentando algo de la propiedad de Insesgadez, sin embargo en comparación con el estimadores 1 y 3, queda en 3er lugar.
Problema 3
Teorema del Límite Central
El Teorema del Límite Central es uno de los más importantes en la inferencia estadística y habla sobre la convergencia de los estimadores como la proporción muestral a la distribución normal. Algunos autores afirman que esta aproximación es bastante buena a partir del umbral \(n > 30\).
A continuación se describen los siguientes pasos para su verificación:
- Realice una simulación en la cual genere una población de \(n=1000\) (Lote), donde el porcentaje de individuos (supongamos plantas) enfermas sea del 50%.
- Genere una función que permita:
- Obtener una muestra aleatoria de la población y
- Calcule el estimador de la proporción muestral \(\widehat{p}\) para un tamaño de muestra dado \(n\).
muestra=function(n){ # esta función saca una muestra de la población
m=sample(poblacion,n,replace=TRUE)
return(m)
}
phat=function(x){ #esta es una función coge la fila de una matriz y saca un promedio de los 1 sobre el tamaño (5), es decir calcula la proporción. Está se integrará dentro de una nueva función para hallar el estimador a partir de varios tamaños de muestra.
sum(x)/5
}- Repita el escenario anterior (b) \(n=500\) veces y analice los resultados en cuanto al comportamiento de los \(500\) resultados del estimador \(\widehat{p}\). ¿Qué tan simétricos o sesgados son los resultados obtenidos? y ¿qué se puede observar en cuanto a la variabilidad?. Realice en su informe un comentario sobre los resultados obtenidos.
Se desarrolla más abajo dentro de los escenarios contemplados en la función, respondiendo así mismo las preguntas.
-
Repita los puntos b y c para tamaños de muestra \(n=5\),
\(10\), \(15\),
\(20\), \(30\),
\(50\), \(60\),
\(100\), \(200\),
\(500\). Compare los resultados
obtenidos para los diferentes tamaños de muestra en cuanto a la
normalidad. Utilice pruebas de bondad y ajuste (shapiro wilks
:
shspiro.test()) y métodos gráficos (gráfico de normalidad:qqnorm()). Comente en su informe los resultados obtenidos
Se crea una función que permita realizar la iteración con cada uno de los escenarios de los tamaños de la muestra, y que adicionalmente permite evaluar el estimador.
grafico1=function(n){ #n es la población.
m=500 # de repeticiones
y=matrix(muestra(n*m),ncol=n) #En muestras de tamaño 5 se calcula la probabilidad
phat=function(x){ #Se ejecuta un estimador de p, esta es una función coje la fila de una matriz y saca un promedio de los 1 sobre el tamaño (5), es decir calcula la proporción
sum(x)/n}
phat5=apply(y,1,phat)
# Creamos el nombre personalizado para el histograma
hist_title <- paste("Histograma para n =", n)
# Creamos el nombre personalizado para el gráfico QQ
qq_title <- paste("Gráfico QQ para n =", n)
# Graficamos el histograma con el nombre personalizado
hist(phat5, main = hist_title)
# Se agrega la visualizacións del gráfico QQ
qqnorm(phat5, main = qq_title)
qqline(phat5)
# Realizamos la prueba de normalidad de Shapiro-Wilk
shapiro_result = shapiro.test(phat5)
# Etiqueta para el resultado del test de Shapiro-Wilk
result_label <- paste("Shapiro-Wilk para n =", n)
print(result_label)
print(shapiro_result)
descr(phat5)
}
par(mfrow=c(2,2))
grafico1(5)## [1] "Shapiro-Wilk para n = 5"
##
## Shapiro-Wilk normality test
##
## data: phat5
## W = 0.93263, p-value = 3.145e-14## Descriptive Statistics
## phat5
## N: 500
##
## phat5
## ----------------- --------
## Mean 0.51
## Std.Dev 0.23
## Min 0.00
## Q1 0.40
## Median 0.60
## Q3 0.60
## Max 1.00
## MAD 0.30
## IQR 0.20
## CV 0.46
## Skewness -0.07
## SE.Skewness 0.11
## Kurtosis -0.43
## N.Valid 500.00
## Pct.Valid 100.00
## [1] "Shapiro-Wilk para n = 10"
##
## Shapiro-Wilk normality test
##
## data: phat5
## W = 0.96267, p-value = 5.913e-10## Descriptive Statistics
## phat5
## N: 500
##
## phat5
## ----------------- --------
## Mean 0.49
## Std.Dev 0.16
## Min 0.00
## Q1 0.40
## Median 0.50
## Q3 0.60
## Max 0.90
## MAD 0.15
## IQR 0.20
## CV 0.32
## Skewness -0.01
## SE.Skewness 0.11
## Kurtosis -0.42
## N.Valid 500.00
## Pct.Valid 100.00## [1] "Shapiro-Wilk para n = 15"
##
## Shapiro-Wilk normality test
##
## data: phat5
## W = 0.97787, p-value = 6.985e-07## Descriptive Statistics
## phat5
## N: 500
##
## phat5
## ----------------- --------
## Mean 0.50
## Std.Dev 0.14
## Min 0.07
## Q1 0.40
## Median 0.53
## Q3 0.60
## Max 0.87
## MAD 0.10
## IQR 0.20
## CV 0.27
## Skewness -0.03
## SE.Skewness 0.11
## Kurtosis -0.24
## N.Valid 500.00
## Pct.Valid 100.00
## [1] "Shapiro-Wilk para n = 20"
##
## Shapiro-Wilk normality test
##
## data: phat5
## W = 0.98117, p-value = 4.549e-06## Descriptive Statistics
## phat5
## N: 500
##
## phat5
## ----------------- --------
## Mean 0.50
## Std.Dev 0.12
## Min 0.10
## Q1 0.40
## Median 0.50
## Q3 0.58
## Max 0.80
## MAD 0.15
## IQR 0.16
## CV 0.23
## Skewness 0.01
## SE.Skewness 0.11
## Kurtosis -0.19
## N.Valid 500.00
## Pct.Valid 100.00## [1] "Shapiro-Wilk para n = 30"
##
## Shapiro-Wilk normality test
##
## data: phat5
## W = 0.98662, p-value = 0.0001497## Descriptive Statistics
## phat5
## N: 500
##
## phat5
## ----------------- --------
## Mean 0.50
## Std.Dev 0.09
## Min 0.20
## Q1 0.43
## Median 0.50
## Q3 0.57
## Max 0.77
## MAD 0.10
## IQR 0.13
## CV 0.19
## Skewness -0.16
## SE.Skewness 0.11
## Kurtosis -0.13
## N.Valid 500.00
## Pct.Valid 100.00
## [1] "Shapiro-Wilk para n = 50"
##
## Shapiro-Wilk normality test
##
## data: phat5
## W = 0.99099, p-value = 0.003848## Descriptive Statistics
## phat5
## N: 500
##
## phat5
## ----------------- --------
## Mean 0.50
## Std.Dev 0.07
## Min 0.26
## Q1 0.44
## Median 0.50
## Q3 0.54
## Max 0.72
## MAD 0.06
## IQR 0.10
## CV 0.15
## Skewness 0.14
## SE.Skewness 0.11
## Kurtosis -0.13
## N.Valid 500.00
## Pct.Valid 100.00## [1] "Shapiro-Wilk para n = 60"
##
## Shapiro-Wilk normality test
##
## data: phat5
## W = 0.99233, p-value = 0.01139## Descriptive Statistics
## phat5
## N: 500
##
## phat5
## ----------------- --------
## Mean 0.50
## Std.Dev 0.06
## Min 0.30
## Q1 0.45
## Median 0.50
## Q3 0.53
## Max 0.70
## MAD 0.05
## IQR 0.08
## CV 0.13
## Skewness -0.03
## SE.Skewness 0.11
## Kurtosis 0.11
## N.Valid 500.00
## Pct.Valid 100.00
## [1] "Shapiro-Wilk para n = 100"
##
## Shapiro-Wilk normality test
##
## data: phat5
## W = 0.99469, p-value = 0.08141## Descriptive Statistics
## phat5
## N: 500
##
## phat5
## ----------------- --------
## Mean 0.50
## Std.Dev 0.05
## Min 0.36
## Q1 0.47
## Median 0.50
## Q3 0.54
## Max 0.65
## MAD 0.04
## IQR 0.07
## CV 0.10
## Skewness 0.02
## SE.Skewness 0.11
## Kurtosis -0.18
## N.Valid 500.00
## Pct.Valid 100.00## [1] "Shapiro-Wilk para n = 200"
##
## Shapiro-Wilk normality test
##
## data: phat5
## W = 0.99645, p-value = 0.3355## Descriptive Statistics
## phat5
## N: 500
##
## phat5
## ----------------- --------
## Mean 0.50
## Std.Dev 0.03
## Min 0.41
## Q1 0.48
## Median 0.50
## Q3 0.52
## Max 0.60
## MAD 0.04
## IQR 0.05
## CV 0.07
## Skewness 0.01
## SE.Skewness 0.11
## Kurtosis -0.18
## N.Valid 500.00
## Pct.Valid 100.00
## [1] "Shapiro-Wilk para n = 500"
##
## Shapiro-Wilk normality test
##
## data: phat5
## W = 0.99757, p-value = 0.6854## Descriptive Statistics
## phat5
## N: 500
##
## phat5
## ----------------- --------
## Mean 0.50
## Std.Dev 0.02
## Min 0.44
## Q1 0.49
## Median 0.50
## Q3 0.52
## Max 0.57
## MAD 0.02
## IQR 0.03
## CV 0.05
## Skewness 0.08
## SE.Skewness 0.11
## Kurtosis -0.03
## N.Valid 500.00
## Pct.Valid 100.00Para determinar la simetría o sesgo de los resultados y observar la variabilidad, podemos analizar varias medidas descriptivas calculadas:
- Se evidencia que a medida que aumenta la muestra, el coeficiente de simetría (Skewness) se acerca a un valor cercano a cero indicando una distribución aproximadamente simétrica, con un sesgo muy pequeños, hacia la derecha o izquierda.
Para analizar la Variabilidad se evidencia:
La desviación estándar (Std.Dev) es decreciente a medida que se aumenta el tamaño de la muestra, lo que indica la dispersión de los datos alrededor de la media, estando cada vez más cercanos a la media.
El coeficiente de variación (CV) tiene a disminuir llegando a valores cercanos a 0, lo que indica una baja variabilidad en relación con la media. Lo mismo sucede con el rango intercuartílico (IQR), lo que indica que la variabilidad en la mitad central de los datos está disminuyendo y que los datos son menos dispersos en torno a la mediana. Esto se debe a una mejor representación de la verdadera variabilidad de la población debido al aumento del tamaño de la muestra.
Test de Shapiro:
En cuanto al test, del estadístico W de manera general podemos analizar cómo inician en un valor 0.93 que asciende hasta 0.99 para la muestra de n=500, en este caso cómo el valor de W se va acercando a 1, sugiere que los datos se ajustan bastante bien a una distribución normal.
Para el valor p, podemos ver valores que inician en p = 3.145e-14 y terminan en p=0.6854 a medida que la muestra incrementa. Esto indica que los datos de la muestra inicial de 5 están lejos de seguir una distribución normal. En este caso, existe evidencia extremadamente fuerte en contra de la hipótesis nula de normalidad.
Por otro lado, el valor p de 0.68 en la prueba de Shapiro-Wilk para la muestra de 500 al final, indicando que no hay suficiente evidencia para rechazar la hipótesis nula de normalidad. En otras palabras, no hay evidencia suficiente para concluir que los datos no provienen de una distribución normal.
La decisión de si considerar los datos como aproximadamente normales o no debe basarse en un análisis integral de los resultados de la prueba y la demás información de contexto disponible.
- Repita toda la simulación (puntos a – d), pero ahora para lotes con 10% de plantas enfermas y de nuevo para lotes con un 90% de plantas enfermas. Concluya sobre los resultados del ejercicio.
Lotes con 10% de plantas enfermas
# Lotes con 10% de plantas enfermas
poblacion10 = c(rep(1,100),rep(0,900))
muestra=function(n){ # esta función saca una muestra de la población
m=sample(poblacion10,n,replace=TRUE)
return(m)
}
grafico=function(n){ #n es la población.
m=500 # de repeticiones
y=matrix(muestra(n*m),ncol=n) #En muestras de tamaño 5 se calcula la probabilidad
phat=function(x){ #Se ejecuta un estimador de p, esta es una función coje la fila de una matriz y saca un promedio de los 1 sobre el tamaño (5), es decir calcula la proporción
sum(x)/n}
phat5=apply(y,1,phat)
# Creamos el nombre personalizado para el histograma
hist_title <- paste("Histograma para n =", n)
# Creamos el nombre personalizado para el gráfico QQ
qq_title <- paste("Gráfico QQ para n =", n)
# Graficamos el histograma con el nombre personalizado
hist(phat5, main = hist_title)
# También podemos agregar otras visualizaciones, como el gráfico QQ
qqnorm(phat5, main = qq_title)
qqline(phat5)
# Realizamos la prueba de normalidad de Shapiro-Wilk
shapiro_result = shapiro.test(phat5)
# Etiqueta para el resultado del test de Shapiro-Wilk
result_label <- paste("Shapiro-Wilk para n =", n)
print(result_label)
print(shapiro_result)
descr(phat5)
}
par(mfrow=c(2,2))
grafico(5)## [1] "Shapiro-Wilk para n = 5"
##
## Shapiro-Wilk normality test
##
## data: phat5
## W = 0.71886, p-value < 2.2e-16## Descriptive Statistics
## phat5
## N: 500
##
## phat5
## ----------------- --------
## Mean 0.10
## Std.Dev 0.13
## Min 0.00
## Q1 0.00
## Median 0.00
## Q3 0.20
## Max 0.60
## MAD 0.00
## IQR 0.20
## CV 1.26
## Skewness 1.02
## SE.Skewness 0.11
## Kurtosis 0.47
## N.Valid 500.00
## Pct.Valid 100.00
## [1] "Shapiro-Wilk para n = 10"
##
## Shapiro-Wilk normality test
##
## data: phat5
## W = 0.84599, p-value < 2.2e-16## Descriptive Statistics
## phat5
## N: 500
##
## phat5
## ----------------- --------
## Mean 0.10
## Std.Dev 0.09
## Min 0.00
## Q1 0.00
## Median 0.10
## Q3 0.20
## Max 0.40
## MAD 0.15
## IQR 0.20
## CV 0.92
## Skewness 0.72
## SE.Skewness 0.11
## Kurtosis -0.05
## N.Valid 500.00
## Pct.Valid 100.00## [1] "Shapiro-Wilk para n = 15"
##
## Shapiro-Wilk normality test
##
## data: phat5
## W = 0.90276, p-value < 2.2e-16## Descriptive Statistics
## phat5
## N: 500
##
## phat5
## ----------------- --------
## Mean 0.10
## Std.Dev 0.08
## Min 0.00
## Q1 0.07
## Median 0.07
## Q3 0.13
## Max 0.40
## MAD 0.10
## IQR 0.07
## CV 0.75
## Skewness 0.54
## SE.Skewness 0.11
## Kurtosis -0.05
## N.Valid 500.00
## Pct.Valid 100.00
## [1] "Shapiro-Wilk para n = 20"
##
## Shapiro-Wilk normality test
##
## data: phat5
## W = 0.91834, p-value = 8.326e-16## Descriptive Statistics
## phat5
## N: 500
##
## phat5
## ----------------- --------
## Mean 0.10
## Std.Dev 0.07
## Min 0.00
## Q1 0.05
## Median 0.10
## Q3 0.15
## Max 0.35
## MAD 0.07
## IQR 0.10
## CV 0.72
## Skewness 0.71
## SE.Skewness 0.11
## Kurtosis 0.43
## N.Valid 500.00
## Pct.Valid 100.00## [1] "Shapiro-Wilk para n = 30"
##
## Shapiro-Wilk normality test
##
## data: phat5
## W = 0.9425, p-value = 5.353e-13## Descriptive Statistics
## phat5
## N: 500
##
## phat5
## ----------------- --------
## Mean 0.10
## Std.Dev 0.05
## Min 0.00
## Q1 0.07
## Median 0.10
## Q3 0.13
## Max 0.33
## MAD 0.05
## IQR 0.07
## CV 0.51
## Skewness 0.62
## SE.Skewness 0.11
## Kurtosis 0.86
## N.Valid 500.00
## Pct.Valid 100.00
## [1] "Shapiro-Wilk para n = 50"
##
## Shapiro-Wilk normality test
##
## data: phat5
## W = 0.96711, p-value = 3.772e-09## Descriptive Statistics
## phat5
## N: 500
##
## phat5
## ----------------- --------
## Mean 0.10
## Std.Dev 0.04
## Min 0.00
## Q1 0.06
## Median 0.10
## Q3 0.12
## Max 0.26
## MAD 0.03
## IQR 0.06
## CV 0.45
## Skewness 0.43
## SE.Skewness 0.11
## Kurtosis 0.04
## N.Valid 500.00
## Pct.Valid 100.00## [1] "Shapiro-Wilk para n = 60"
##
## Shapiro-Wilk normality test
##
## data: phat5
## W = 0.98149, p-value = 5.523e-06## Descriptive Statistics
## phat5
## N: 500
##
## phat5
## ----------------- --------
## Mean 0.10
## Std.Dev 0.04
## Min 0.00
## Q1 0.07
## Median 0.10
## Q3 0.12
## Max 0.23
## MAD 0.05
## IQR 0.05
## CV 0.40
## Skewness 0.20
## SE.Skewness 0.11
## Kurtosis 0.14
## N.Valid 500.00
## Pct.Valid 100.00
## [1] "Shapiro-Wilk para n = 100"
##
## Shapiro-Wilk normality test
##
## data: phat5
## W = 0.98143, p-value = 5.328e-06## Descriptive Statistics
## phat5
## N: 500
##
## phat5
## ----------------- --------
## Mean 0.10
## Std.Dev 0.03
## Min 0.02
## Q1 0.08
## Median 0.09
## Q3 0.12
## Max 0.19
## MAD 0.03
## IQR 0.04
## CV 0.30
## Skewness 0.32
## SE.Skewness 0.11
## Kurtosis -0.04
## N.Valid 500.00
## Pct.Valid 100.00## [1] "Shapiro-Wilk para n = 200"
##
## Shapiro-Wilk normality test
##
## data: phat5
## W = 0.99166, p-value = 0.006594## Descriptive Statistics
## phat5
## N: 500
##
## phat5
## ----------------- --------
## Mean 0.10
## Std.Dev 0.02
## Min 0.04
## Q1 0.09
## Median 0.10
## Q3 0.12
## Max 0.17
## MAD 0.02
## IQR 0.03
## CV 0.21
## Skewness 0.11
## SE.Skewness 0.11
## Kurtosis -0.03
## N.Valid 500.00
## Pct.Valid 100.00
## [1] "Shapiro-Wilk para n = 500"
##
## Shapiro-Wilk normality test
##
## data: phat5
## W = 0.9941, p-value = 0.04982## Descriptive Statistics
## phat5
## N: 500
##
## phat5
## ----------------- --------
## Mean 0.10
## Std.Dev 0.01
## Min 0.06
## Q1 0.09
## Median 0.10
## Q3 0.11
## Max 0.14
## MAD 0.01
## IQR 0.02
## CV 0.13
## Skewness 0.16
## SE.Skewness 0.11
## Kurtosis -0.23
## N.Valid 500.00
## Pct.Valid 100.00Lotes con 90% de plantas enfermas
# Lotes con 90% de plantas enfermas
poblacion90 = c(rep(1,900),rep(0,100))
muestra=function(n){ # esta función saca una muestra de la población
m=sample(poblacion90,n,replace=TRUE)
return(m)
}
grafico3=function(n){ #n es la población.
m=500 # de repeticiones
y=matrix(muestra(n*m),ncol=n) #En muestras de tamaño 5 se calcula la probabilidad
phat=function(x){ #Se ejecuta un estimador de p, esta es una función coje la fila de una matriz y saca un promedio de los 1 sobre el tamaño (5), es decir calcula la proporción
sum(x)/n}
phat5=apply(y,1,phat)
# Creamos el nombre personalizado para el histograma
hist_title <- paste("Histograma para n =", n)
# Creamos el nombre personalizado para el gráfico QQ
qq_title <- paste("Gráfico QQ para n =", n)
# Graficamos el histograma con el nombre personalizado
hist(phat5, main = hist_title)
# También podemos agregar otras visualizaciones, como el gráfico QQ
qqnorm(phat5, main = qq_title)
qqline(phat5)
# Realizamos la prueba de normalidad de Shapiro-Wilk
shapiro_result = shapiro.test(phat5)
# Etiqueta para el resultado del test de Shapiro-Wilk
result_label <- paste("Shapiro-Wilk para n =", n)
print(result_label)
print(shapiro_result)
descr(phat5)
}
par(mfrow=c(2,2))
grafico3(5)## [1] "Shapiro-Wilk para n = 5"
##
## Shapiro-Wilk normality test
##
## data: phat5
## W = 0.69394, p-value < 2.2e-16## Descriptive Statistics
## phat5
## N: 500
##
## phat5
## ----------------- --------
## Mean 0.91
## Std.Dev 0.12
## Min 0.40
## Q1 0.80
## Median 1.00
## Q3 1.00
## Max 1.00
## MAD 0.00
## IQR 0.20
## CV 0.14
## Skewness -1.10
## SE.Skewness 0.11
## Kurtosis 0.62
## N.Valid 500.00
## Pct.Valid 100.00
## [1] "Shapiro-Wilk para n = 10"
##
## Shapiro-Wilk normality test
##
## data: phat5
## W = 0.84263, p-value < 2.2e-16## Descriptive Statistics
## phat5
## N: 500
##
## phat5
## ----------------- --------
## Mean 0.90
## Std.Dev 0.10
## Min 0.60
## Q1 0.80
## Median 0.90
## Q3 1.00
## Max 1.00
## MAD 0.15
## IQR 0.20
## CV 0.11
## Skewness -0.74
## SE.Skewness 0.11
## Kurtosis -0.18
## N.Valid 500.00
## Pct.Valid 100.00## [1] "Shapiro-Wilk para n = 15"
##
## Shapiro-Wilk normality test
##
## data: phat5
## W = 0.8839, p-value < 2.2e-16## Descriptive Statistics
## phat5
## N: 500
##
## phat5
## ----------------- --------
## Mean 0.91
## Std.Dev 0.08
## Min 0.67
## Q1 0.87
## Median 0.93
## Q3 1.00
## Max 1.00
## MAD 0.10
## IQR 0.13
## CV 0.09
## Skewness -0.68
## SE.Skewness 0.11
## Kurtosis -0.24
## N.Valid 500.00
## Pct.Valid 100.00
## [1] "Shapiro-Wilk para n = 20"
##
## Shapiro-Wilk normality test
##
## data: phat5
## W = 0.92118, p-value = 1.649e-15## Descriptive Statistics
## phat5
## N: 500
##
## phat5
## ----------------- --------
## Mean 0.89
## Std.Dev 0.07
## Min 0.65
## Q1 0.85
## Median 0.90
## Q3 0.95
## Max 1.00
## MAD 0.07
## IQR 0.10
## CV 0.08
## Skewness -0.70
## SE.Skewness 0.11
## Kurtosis 0.51
## N.Valid 500.00
## Pct.Valid 100.00## [1] "Shapiro-Wilk para n = 30"
##
## Shapiro-Wilk normality test
##
## data: phat5
## W = 0.95019, p-value = 6.184e-12## Descriptive Statistics
## phat5
## N: 500
##
## phat5
## ----------------- --------
## Mean 0.90
## Std.Dev 0.06
## Min 0.70
## Q1 0.87
## Median 0.90
## Q3 0.93
## Max 1.00
## MAD 0.05
## IQR 0.07
## CV 0.06
## Skewness -0.49
## SE.Skewness 0.11
## Kurtosis 0.21
## N.Valid 500.00
## Pct.Valid 100.00
## [1] "Shapiro-Wilk para n = 50"
##
## Shapiro-Wilk normality test
##
## data: phat5
## W = 0.96889, p-value = 8.294e-09## Descriptive Statistics
## phat5
## N: 500
##
## phat5
## ----------------- --------
## Mean 0.90
## Std.Dev 0.04
## Min 0.76
## Q1 0.88
## Median 0.90
## Q3 0.92
## Max 1.00
## MAD 0.03
## IQR 0.04
## CV 0.05
## Skewness -0.38
## SE.Skewness 0.11
## Kurtosis -0.10
## N.Valid 500.00
## Pct.Valid 100.00## [1] "Shapiro-Wilk para n = 60"
##
## Shapiro-Wilk normality test
##
## data: phat5
## W = 0.97176, p-value = 3.118e-08## Descriptive Statistics
## phat5
## N: 500
##
## phat5
## ----------------- --------
## Mean 0.90
## Std.Dev 0.04
## Min 0.77
## Q1 0.88
## Median 0.90
## Q3 0.93
## Max 1.00
## MAD 0.05
## IQR 0.05
## CV 0.04
## Skewness -0.42
## SE.Skewness 0.11
## Kurtosis 0.20
## N.Valid 500.00
## Pct.Valid 100.00
## [1] "Shapiro-Wilk para n = 100"
##
## Shapiro-Wilk normality test
##
## data: phat5
## W = 0.98289, p-value = 1.294e-05## Descriptive Statistics
## phat5
## N: 500
##
## phat5
## ----------------- --------
## Mean 0.90
## Std.Dev 0.03
## Min 0.80
## Q1 0.88
## Median 0.90
## Q3 0.92
## Max 0.98
## MAD 0.03
## IQR 0.04
## CV 0.03
## Skewness -0.30
## SE.Skewness 0.11
## Kurtosis -0.01
## N.Valid 500.00
## Pct.Valid 100.00## [1] "Shapiro-Wilk para n = 200"
##
## Shapiro-Wilk normality test
##
## data: phat5
## W = 0.98912, p-value = 0.0009086## Descriptive Statistics
## phat5
## N: 500
##
## phat5
## ----------------- --------
## Mean 0.90
## Std.Dev 0.02
## Min 0.83
## Q1 0.88
## Median 0.90
## Q3 0.92
## Max 0.96
## MAD 0.02
## IQR 0.03
## CV 0.02
## Skewness -0.23
## SE.Skewness 0.11
## Kurtosis 0.35
## N.Valid 500.00
## Pct.Valid 100.00
## [1] "Shapiro-Wilk para n = 500"
##
## Shapiro-Wilk normality test
##
## data: phat5
## W = 0.99231, p-value = 0.0112## Descriptive Statistics
## phat5
## N: 500
##
## phat5
## ----------------- --------
## Mean 0.90
## Std.Dev 0.01
## Min 0.85
## Q1 0.89
## Median 0.90
## Q3 0.91
## Max 0.93
## MAD 0.01
## IQR 0.02
## CV 0.02
## Skewness -0.08
## SE.Skewness 0.11
## Kurtosis -0.29
## N.Valid 500.00
## Pct.Valid 100.00Realizando el ejercicio ejercicio pero condiciones particulares en las muestras diferentes, se ve cómo las estadísticas de las muestras, como la media, la mediana, la varianza y la asimetría, cambian a medida que se aumenta el tamaño de la muestra, y estos cambios van en función con lo esperado en el Teorema del Límite Central demostrandose en su totalidad.
Un punto relvante a mencionar es que si analizamos las gráficas de histograma y el qq-normalidad, la diferencia entre los dos casos es que inician con una distribución sesgada hacia un lado determinado y con vacios visualmente apreciables al partir de datos discretos y muestras pequeñas. En el caso de 10% Plantas enfermas, la distribución es identifica sesgada hacia la izquierda, mientras que con el 90% Plantas enfermas, la distribución es sesgada hacia la derecha.
En los tres casos analizados, cualquiera que sea la distribución de la proporción de plantas enfermas de la población, a medida en que aumentan el tamaño de la muestra, se distribuye de manera normal y la estimación tambien se acerca con mayor precisión al valor real,por lo tanto se confirma el teorema del límite central.
Problema 4
Estimacción boostrap
Cuando se extrae una muestra de una población que no es normal y se requiere estimar un intervalo de confianza se pueden utilizar los métodos de estimación bootstrap. Esta metodología supone que se puede reconstruir la población objeto de estudio mediante un muestreo con reemplazo de la muestra que se tiene. Existen varias versiones del método. Una presentación básica del método se describe a continuación:
El artículo de In-use Emissions from Heavy Duty Dissel Vehicles (J.Yanowitz, 2001) presenta las mediciones de eficiencia de combustible en millas/galón de una muestra de siete camiones. Los datos obtenidos son los siguientes: 7.69, 4.97, 4.56, 6.49, 4.34, 6.24 y 4.45. Se supone que es una muestra aleatoria de camiones y que se desea construir un intervalo de confianza del 95 % para la media de la eficiencia de combustible de esta población. No se tiene información de la distribución de los datos. El método bootstrap permite construir intervalos de confianza del 95 % - Para ilustrar el método suponga que coloca los valores de la muestra en una caja y extrae uno al azar. Este correspondería al primer valor de la muestra bootstrap \(X^{∗}_{1}\). Después de anotado el valor se regresa \(X^{∗}_{1}\) a la caja y se extrae el valor \(X^{∗}_{2}\) , regresandolo nuevamente. Este procedimiento se repite hasta completar una muestra de tamaño \(n\), \(X^{∗}_{1}\),\(X^{∗}_{2}\),\(X^{∗}_{2}\),\(X^{∗}_{n}\), conformando la muestra bootstrap.
Es necesario extraer un gran número de muestras (suponga k = 1000). Para cada una de las muestra bootstrap obtenidas se calcula la media \(\bar{X^{∗}_{i}}\), obteniéndose un valor para cada muestra. El intervalo de confianza queda conformado por los percentiles \(P_{2.5}\) y \(P_{97.5}\). Existen dos métodos para estimarlo:
| Método 1 | \((P_{2.5}; P_{97.5})\) |
| Método 2 | \((2\bar{X} − P_{97.5}; 2\bar{X} − P_{2.5})\) |
Construya el intervalo de confianza por los dos métodos y compare los resultados obtenidos. Comente los resultados. Confiaría en estas estimaciones?
set.seed(123)
muestra=c( 7.69, 4.97, 4.56, 6.49, 4.34, 6.24, 4.45) # datos muestra
boot=sample(muestra,7000,replace=TRUE) # se extraen n x m muestras
mtrx=matrix(boot,nrow=1000,ncol=7) # se construye matriz de n x m
md=apply(mtrx,1,mean) # se calculan m medias por cada fila de esa matriz, md es un valor de 1000 medias
# Forma 1
ic1=quantile(md, probs=c(0.025, 0.975)) # A las 1000 medias se le saca los percentiles 25 y 97.5. buscando tener el intervalo de confianza del 95%, se calcula IC método 1
ic1## 2.5% 97.5%
## 4.721429 6.418643# Forma 2
ic2=c(2*mean(md)-ic1[2], 2*mean(md)-ic1[1]) # Se toma el vector, la media restando el límite superior y 2 veces la media menos el límite inferior. se calcula IC método 2
ic2## 97.5% 2.5%
## 4.639526 6.336740colores <- c("#519DE9", "#4CB140", "#C9190B")
# Se crea histograma para representar la estimación boostrap.
hist(md, las = 1, main = "Histograma", ylab = "Frecuencia", xlab = "Valores", col = colores[1])
# Se trazan unas líneas verticales para los intervalos de confianza.
abline(v = ic1, col = colores[2], lwd = 2)
abline(v = ic2, col = colores[3], lwd = 2)
- En este caso se está buscando estimar la eficiencia de combustible en millas por galón a partir del muestreo de 7 camiones.
- Interpretando los resultados con el método de estimación no parámetrica # 1, el promedio de eficiencia de gasolina por millas por galón se encuentra entre 4.72 y 6.42 con un 95% de confianza. Mientras que en el 2do método, está entre 4.64 y 6.34 millas por galón bajo el mismo intervalo de confianza.