Actividad 2

Problema 1: Estimación del valor de π

En la figura, un círculo con un área igual a π/4 está inscrito en un cuadrado cuya área es igual a 1. Se eligen 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 éste, la cual es π/4, por tanto, se puede estimar el valor de π/4 al contar el número de puntos dentro del círculo, para obtener así una aproximación para el valor de π.

¿Cuántos de los puntos están dentro del círculo? ¿Cuál es su estimación de π?

#

n <- 100000
x <- runif(n, min = 0, max = 1)
y <- runif(n, min = 0, max = 1)
dc <- (x - 0.5)^2 + (y - 0.5)^2 #Distancia desde el centro
eval <- as.numeric(dc <= 0.25) #Evaluar si dc <= 0.5
coordinates <- data.frame(x, y, dc, eval)
dentro <- marginSums(coordinates$eval) #Obtener puntos dc <= 0.5
prob <- dentro/n #Prob de tener un punto dc <= 0.5
estpi <- prob*4 #Estimacion de pi
error <- pi - estpi
cat("Pares de coordenadas:", n,
"\nPuntos con distancia al centro menor a 0.5:", dentro,
"\nProporción de puntos dentro del círculo:", prob,
"\nValor estimado de pi:", estpi,
"\nError:", error)

Pares de coordenadas: 1e+05 
Puntos con distancia al centro menor a 0.5: 78792 
Proporción de puntos dentro del círculo: 0.78792 
Valor estimado de pi: 3.15168 
Error: -0.01008735

Problema 2: Propiedades de los estimadores

Sean X1 ,X2 ,X3 y X4, una muestra aleatoria de tamaño n=4 cuya población la conforma una distribución exponencial con parámetro θ desconocido. Determine las características de cada uno de los siguientes estimadores propuestos:

#

Genere una muestras de n=20, 50, 100 y 1000 para cada uno de los estimadores planteados.

En cada caso evalue las propiedades de insesgadez, eficiencia y consistencia. Suponga un valor para el parámetro θ

Tamaño de muestra n=20

n <- 20
param <- 5
x1 <- rexp(n, rate = 1/param)
x2 <- rexp(n, rate = 1/param)
x3 <- rexp(n, rate = 1/param)
x4 <- rexp(n, rate = 1/param)

estm1 <- (x1 + x2)/6 + (x3 + x4)/3
estm2 <- (x1 + 2*x2 + 3*x3 + 4*x4)/5
estm3 <- (x1 + x2 + x3 + x4)/4
estm4 <- 0
for (i in 1:n) {
  estm4[i] <- (min(x1[i], x2[i], x3[i], x4[i]) + max(x1[i], x2[i], x3[i], x4[i]))/2
}

data <- data.frame(estm1, estm2, estm3, estm4)
boxplot(data, las = 1)
abline(h = param, col = "red")

media <- apply(data, 2, mean)
varianza <- apply(data, 2, var)

rbind(media, varianza)

            estm1    estm2     estm3     estm4
media     6.02333 12.36359  5.712452  6.896541
varianza 14.01798 64.18858 12.498340 18.358600

Se observa que los mejores resultados se obtienen para los estimadores 1, 3 y 4, ya que tienen una media más cercana al parámetro. De manera contraria, se observa que el estimador 2 es un estimador sesgado. De los estimadores que presentan mejor comportamiento, el estm3 tiene la menor desviación estándar, por lo que sería el estimador eficiente.

Tamaño de muestra n=50

n <- 50
param <- 5
x1 <- rexp(n, rate = 1/param)
x2 <- rexp(n, rate = 1/param)
x3 <- rexp(n, rate = 1/param)
x4 <- rexp(n, rate = 1/param)

estm1 <- (x1 + x2)/6 + (x3 + x4)/3
estm2 <- (x1 + 2*x2 + 3*x3 + 4*x4)/5
estm3 <- (x1 + x2 + x3 + x4)/4
estm4 <- 0
for (i in 1:n) {
  estm4[i] <- (min(x1[i], x2[i], x3[i], x4[i]) + max(x1[i], x2[i], x3[i], x4[i]))/2
}

data <- data.frame(estm1, estm2, estm3, estm4)
boxplot(data, las = 1)
abline(h = param, col = "red")

media <- apply(data, 2, mean)
varianza <- apply(data, 2, var)

rbind(media, varianza)

            estm1    estm2    estm3    estm4
media    5.424643 11.08826 5.374657 6.271237
varianza 5.655893 23.10447 4.924779 8.335489

Tamaño de muestra n=100

n <- 100
param <- 5
x1 <- rexp(n, rate = 1/param)
x2 <- rexp(n, rate = 1/param)
x3 <- rexp(n, rate = 1/param)
x4 <- rexp(n, rate = 1/param)

estm1 <- (x1 + x2)/6 + (x3 + x4)/3
estm2 <- (x1 + 2*x2 + 3*x3 + 4*x4)/5
estm3 <- (x1 + x2 + x3 + x4)/4
estm4 <- 0
for (i in 1:n) {
  estm4[i] <- (min(x1[i], x2[i], x3[i], x4[i]) + max(x1[i], x2[i], x3[i], x4[i]))/2
}

data <- data.frame(estm1, estm2, estm3, estm4)
boxplot(data, las = 1)
abline(h = param, col = "red")

media <- apply(data, 2, mean)
varianza <- apply(data, 2, var)

rbind(media, varianza)

            estm1    estm2    estm3     estm4
media    5.365317 10.62637 5.334601  6.108345
varianza 7.946827 32.69939 7.253206 11.365050

Tamaño de muestra n=1000

n <- 1000
param <- 5
x1 <- rexp(n, rate = 1/param)
x2 <- rexp(n, rate = 1/param)
x3 <- rexp(n, rate = 1/param)
x4 <- rexp(n, rate = 1/param)

estm1 <- (x1 + x2)/6 + (x3 + x4)/3
estm2 <- (x1 + 2*x2 + 3*x3 + 4*x4)/5
estm3 <- (x1 + x2 + x3 + x4)/4
estm4 <- 0
for (i in 1:n) {
  estm4[i] <- (min(x1[i], x2[i], x3[i], x4[i]) + max(x1[i], x2[i], x3[i], x4[i]))/2
}

data <- data.frame(estm1, estm2, estm3, estm4)
boxplot(data, las = 1)
abline(h = param, col = "red")

media <- apply(data, 2, mean)
varianza <- apply(data, 2, var)

rbind(media, varianza)

            estm1     estm2    estm3     estm4
media    4.874131  9.759504 4.887361  5.666101
varianza 7.025220 30.163163 6.371407 10.739396

A medida que aumenta el número de muestras, se observa que los mejores resultados continuan siendo para los estimadores 1, 3 y 4, mientras que el sesgo del 2 disminuye. El estimador 3 continúa siendo el eficiente. Finalmente, se observa que al aumentar la muestra, la varianza de los indicadores se reduce.

Problema 3: Teorema de límite central

Realice una simulación en la cual genere una población N=1000, donde el porcentaje de individuos (supongamos plantas) enfermas sea del 50%.

# Generación de la población simulada
N <- 1000 # Tamaño poblacion
P <- 0.5 # parámetro proporción de plantas enfermas
poblacion <- rbinom(N, 1, P) # Población de plantas

Genere una función que permita obtener una muestra aleatoria de la población y calcule el estimador de la proporción muestral p para un tamaño de muestra dado n.

# Generación del estimador de la proporción
n <- 500 # Tamaño muestra
muestra <- sample(poblacion, n) # muestra de n plantas
tablmuestra <- table(muestra)
propmuestra <- prop.table(tablmuestra)
p <- sum(muestra) / n # Estimador de la proporción

Repita el escenario anterior 500 veces y analice los resultados en cuanto al comportamiento de los resultados del estimador para tamaños de muestra n=5, n=10, n=15, n=20, n=30, n=50, n=60, n=100, n=200 y n=500 comparando 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

¿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.

Análisis de estimadores variando la muestra

# Generación del estimador 500 veces
p <- 0
Estimar500 <- function(n) {
    for (i in 1:500) {
        muestra <- sample(poblacion, n)
        p[i] <- sum(muestra) / n
    }
    return(p)
}

p5 <- Estimar500(5)
p10 <- Estimar500(10)
p15 <- Estimar500(15)
p20 <- Estimar500(20)
p30 <- Estimar500(30)
p50 <- Estimar500(50)
p60 <- Estimar500(60)
p100 <- Estimar500(100)
p200 <- Estimar500(200)
p500 <- Estimar500(500)

par(mfrow = c(2, 5))
boxplot(p5, xlab = "n = 5", las = 2)
abline(h = 0.5, col = "red")
boxplot(p10, xlab = "n = 10", las = 2)
abline(h = 0.5, col = "red")
boxplot(p15, xlab = "n = 15", las = 2)
abline(h = 0.5, col = "red")
boxplot(p20, xlab = "n = 20", las = 2)
abline(h = 0.5, col = "red")
boxplot(p30, xlab = "n = 30", las = 2)
abline(h = 0.5, col = "red")
boxplot(p50, xlab = "n = 50", las = 2)
abline(h = 0.5, col = "red")
boxplot(p60, xlab = "n = 60", las = 2)
abline(h = 0.5, col = "red")
boxplot(p100, xlab = "n = 100", las = 2)
abline(h = 0.5, col = "red")
boxplot(p200, xlab = "n = 200", las = 2)
abline(h = 0.5, col = "red")
boxplot(p500, xlab = "n = 500", las = 2)
abline(h = 0.5, col = "red")

analisisp5 <- rbind(round(mean(p5),3), round(var(p5),3))
analisisp10 <- rbind(round(mean(p10),3), round(var(p10),3))
analisisp15 <- rbind(round(mean(p15),3), round(var(p15),3))
analisisp20 <- rbind(round(mean(p20),3), round(var(p20),3))
analisisp30 <- rbind(round(mean(p30),3), round(var(p30),3))
analisisp50 <- rbind(round(mean(p50),3), round(var(p50),3))
analisisp60 <- rbind(round(mean(p60),3), round(var(p60),3))
analisisp100 <- rbind(round(mean(p100),3), round(var(p100),3))
analisisp200 <- rbind(round(mean(p200),3), round(var(p200),3))
analisisp500 <- rbind(round(mean(p500),3), round(var(p500),3))
analsisisestimadores <- cbind(analisisp5, analisisp10, analisisp15, analisisp20, analisisp30, analisisp50,
analisisp60, analisisp100, analisisp200, analisisp500)
colnames(analsisisestimadores) <- c("n=5", "n=10", "n=15", "n=20", "n=30", "n=50", "n=60", "n=100", "n=200", "n=500")
rownames(analsisisestimadores) <- c("media", "varianza")
print(analsisisestimadores)

           n=5  n=10  n=15  n=20  n=30  n=50  n=60 n=100 n=200 n=500
media    0.502 0.495 0.513 0.501 0.506 0.508 0.503 0.504 0.504 0.506
varianza 0.056 0.024 0.014 0.012 0.009 0.005 0.004 0.002 0.001 0.000

Se observa que el estimador que presenta el menor sesgo, es el calculado con una muestra n=20 repetida 500 veces, de igual manera, es el estimador cuyo diagrama de caja presenta menor longitud, por lo que vendría a ser un estimador eficiente. Esto demuestra que la calidad de los estimadores puede ser muy buena sin necesidad de tener un número de muestras excesivo, por lo que, en la planeación de los proyectos es muy importante determinar un tamaño de muestra adecuado con el fin de no perder eficiencia sin la necesidad de realizar un número de muestras excesivo que pueda aumentar el costo del mismo.

Análisis de normalidad variando la muestra

En el siguiente arreglo de gráficos, se presenta un histograma del estimador para varios tamaños de muestra.

# Generación del estimador 500 veces
par(mfrow = c(2, 5))
hist(p5)
abline(v = mean(p5), col = "red")
abline(v = median(p5), col = "blue")
hist(p10)
abline(v = mean(p10), col = "red")
abline(v = median(p10), col = "blue")
hist(p15)
abline(v = mean(p15), col = "red")
abline(v = median(p15), col = "blue")
hist(p20)
abline(v = mean(p20), col = "red")
abline(v = median(p20), col = "blue")
hist(p30)
abline(v = mean(p30), col = "red")
abline(v = median(p30), col = "blue")
hist(p50)
abline(v = mean(p50), col = "red")
abline(v = median(p50), col = "blue")
hist(p60)
abline(v = mean(p60), col = "red")
abline(v = median(p60), col = "blue")
hist(p100)
abline(v = mean(p100), col = "red")
abline(v = median(p100), col = "blue")
hist(p200)
abline(v = mean(p200), col = "red")
abline(v = median(p200), col = "blue")
hist(p500)
abline(v = mean(p500), col = "red")
abline(v = median(p500), col = "blue")

Lo que se observa es que al aumentar el tamaño de muestra, la distribución tiende a comportarse de manera simétrica, pues la media y la moda (linea roja y azul respectivamente) se van acercando cada vez más.

Otra método gráfico para observar este teorema, son los gráficos qq presentados a continuación.

par(mfrow = c(2, 5))
qqnorm(p5, main = "Normal Q-Q n=5")
qqline(p5, col = "red")
qqnorm(p10, main = "Normal Q-Q n=10")
qqline(p10, col = "red")
qqnorm(p15, main = "Normal Q-Q n=15")
qqline(p15, col = "red")
qqnorm(p20, main = "Normal Q-Q n=20")
qqline(p20, col = "red")
qqnorm(p30, main = "Normal Q-Q n=30")
qqline(p30, col = "red")
qqnorm(p50, main = "Normal Q-Q  n=50")
qqline(p50, col = "red")
qqnorm(p60, main = "Normal Q-Q n=60")
qqline(p60, col = "red")
qqnorm(p100, main = "Normal Q-Q n=100")
qqline(p100, col = "red")
qqnorm(p200, main = "Normal Q-Q n=200")
qqline(p200, col = "red")
qqnorm(p500, main = "Normal Q-Q n=500")
qqline(p500, col = "red")

Igualmente, al aumentar el tamaño de la muestra, la distribución tiende a comportarse como una distribución normal, pues los quantiles reales tienden a tomar el valor teórico de una distribución normal.

Finalmente, se contrastaron dos tests de normalidad con una Ho: distribución normal y Ha: no distribución normal. Los resultados se presentan a continuación.

library(nortest)
library(formattable)
options(digits = 2)
shapiro <- c(shapiro.test(p5)$p, shapiro.test(p10)$p, shapiro.test(p15)$p,
              shapiro.test(p20)$p, shapiro.test(p30)$p, shapiro.test(p50)$p,
              shapiro.test(p60)$p, shapiro.test(p100)$p, shapiro.test(p200)$p,
              shapiro.test(p500)$p)
lillie <- c(lillie.test(p5)$p, lillie.test(p10)$p, lillie.test(p15)$p,
              lillie.test(p20)$p, lillie.test(p30)$p, lillie.test(p50)$p,
              lillie.test(p60)$p, lillie.test(p100)$p, lillie.test(p200)$p,
              lillie.test(p500)$p)
pvalue <- data.frame(shapiro, lillie)
rownames(pvalue) <- c("n5", "n10", "n15", "n20", "n30", "n50", "n60",
                        "n100", "n200", "n500")

library(formattable)
formattable(pvalue, list(
  shapiro = formatter("span",
    style = x ~ style(color = ifelse(pvalue$shapiro < 0.05, "red", "green"))),
  lillie = formatter("span",
    style = x ~ style(color = ifelse(pvalue$lillie < 0.05, "red", "green")))))

	shapiro	lillie
n5	3.2e-14	8.7e-47
n10	1.0e-09	4.7e-23
n15	1.6e-08	4.4e-18
n20	4.9e-07	6.5e-15
n30	2.2e-04	1.7e-07
n50	1.0e-02	1.9e-05
n60	3.9e-03	1.6e-05
n100	1.0e-01	4.8e-03
n200	8.3e-02	5.1e-03
n500	8.9e-02	1.1e-01

Se observa que con la prueba de Shapiro-Wilk no se rechaza la hipótesis nula a partir de la simulación realizada con la toma de 60 muestras en 500 ocasiones, contrariamente, con la prueba de Kolmogorov-Smirnov (Lilliefors) no se rechaza en la última simulación, cuando se toman 500 muestras para calcular la proporción muestral en 500 ocasiones.

La razón por la que esto podría suceder, es que la prueba de Shapiro-Wilk con la variación de Patrick Royston en 1982 adecua el uso del método para muestras entre 3 y 2000 datos (Royston, J. P., 1982), mientras que la prueba de Lilliefors se recomienda para un número de muestras mayores a 50. En el caso de este ejercicio, se calcularon 500 valores diferentes del estimador de la proporción, por lo que los reusltados de Shapiro-Wilk con la variación mencionada son menos conservadores que la de Lillifors.

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.

10% Plantas enfermas

P <- 0.1 # parámetro proporción de plantas enfermas
poblacion <- rbinom(N, 1, P) # Población de plantas

p5 <- Estimar500(5)
p10 <- Estimar500(10)
p15 <- Estimar500(15)
p20 <- Estimar500(20)
p30 <- Estimar500(30)
p50 <- Estimar500(50)
p60 <- Estimar500(60)
p100 <- Estimar500(100)
p200 <- Estimar500(200)
p500 <- Estimar500(500)

par(mfrow = c(2, 5))
boxplot(p5, xlab = "n = 5", las = 2)
abline(h = 0.1, col = "red")
boxplot(p10, xlab = "n = 10", las = 2)
abline(h = 0.1, col = "red")
boxplot(p15, xlab = "n = 15", las = 2)
abline(h = 0.1, col = "red")
boxplot(p20, xlab = "n = 20", las = 2)
abline(h = 0.1, col = "red")
boxplot(p30, xlab = "n = 30", las = 2)
abline(h = 0.1, col = "red")
boxplot(p50, xlab = "n = 50", las = 2)
abline(h = 0.1, col = "red")
boxplot(p60, xlab = "n = 60", las = 2)
abline(h = 0.1, col = "red")
boxplot(p100, xlab = "n = 100", las = 2)
abline(h = 0.1, col = "red")
boxplot(p200, xlab = "n = 200", las = 2)
abline(h = 0.1, col = "red")
boxplot(p500, xlab = "n = 500", las = 2)
abline(h = 0.1, col = "red")

analisisp5 <- rbind(round(mean(p5),3), round(var(p5),3))
analisisp10 <- rbind(round(mean(p10),3), round(var(p10),3))
analisisp15 <- rbind(round(mean(p15),3), round(var(p15),3))
analisisp20 <- rbind(round(mean(p20),3), round(var(p20),3))
analisisp30 <- rbind(round(mean(p30),3), round(var(p30),3))
analisisp50 <- rbind(round(mean(p50),3), round(var(p50),3))
analisisp60 <- rbind(round(mean(p60),3), round(var(p60),3))
analisisp100 <- rbind(round(mean(p100),3), round(var(p100),3))
analisisp200 <- rbind(round(mean(p200),3), round(var(p200),3))
analisisp500 <- rbind(round(mean(p500),3), round(var(p500),3))
analsisisestimadores <- cbind(analisisp5, analisisp10, analisisp15, analisisp20, analisisp30, analisisp50,
analisisp60, analisisp100, analisisp200, analisisp500)
colnames(analsisisestimadores) <- c("n=5", "n=10", "n=15", "n=20", "n=30", "n=50", "n=60", "n=100", "n=200", "n=500")
rownames(analsisisestimadores) <- c("media", "varianza")
print(analsisisestimadores)

           n=5  n=10  n=15  n=20  n=30  n=50  n=60 n=100 n=200 n=500
media    0.123 0.123 0.119 0.124 0.123 0.123 0.124 0.124  0.12  0.12
varianza 0.022 0.011 0.007 0.005 0.003 0.002 0.002 0.001  0.00  0.00

par(mfrow = c(2, 5))
hist(p5)
abline(v = mean(p5), col = "red")
abline(v = median(p5), col = "blue")
hist(p10)
abline(v = mean(p10), col = "red")
abline(v = median(p10), col = "blue")
hist(p15)
abline(v = mean(p15), col = "red")
abline(v = median(p15), col = "blue")
hist(p20)
abline(v = mean(p20), col = "red")
abline(v = median(p20), col = "blue")
hist(p30)
abline(v = mean(p30), col = "red")
abline(v = median(p30), col = "blue")
hist(p50)
abline(v = mean(p50), col = "red")
abline(v = median(p50), col = "blue")
hist(p60)
abline(v = mean(p60), col = "red")
abline(v = median(p60), col = "blue")
hist(p100)
abline(v = mean(p100), col = "red")
abline(v = median(p100), col = "blue")
hist(p200)
abline(v = mean(p200), col = "red")
abline(v = median(p200), col = "blue")
hist(p500)
abline(v = mean(p500), col = "red")
abline(v = median(p500), col = "blue")

par(mfrow = c(2, 5))
qqnorm(p5, main = "Normal Q-Q n=5")
qqline(p5, col = "red")
qqnorm(p10, main = "Normal Q-Q n=10")
qqline(p10, col = "red")
qqnorm(p15, main = "Normal Q-Q n=15")
qqline(p15, col = "red")
qqnorm(p20, main = "Normal Q-Q n=20")
qqline(p20, col = "red")
qqnorm(p30, main = "Normal Q-Q n=30")
qqline(p30, col = "red")
qqnorm(p50, main = "Normal Q-Q  n=50")
qqline(p50, col = "red")
qqnorm(p60, main = "Normal Q-Q n=60")
qqline(p60, col = "red")
qqnorm(p100, main = "Normal Q-Q n=100")
qqline(p100, col = "red")
qqnorm(p200, main = "Normal Q-Q n=200")
qqline(p200, col = "red")
qqnorm(p500, main = "Normal Q-Q n=500")
qqline(p500, col = "red")

library(nortest)
library(formattable)
options(digits = 2)
shapiro <- c(shapiro.test(p5)$p, shapiro.test(p10)$p, shapiro.test(p15)$p,
              shapiro.test(p20)$p, shapiro.test(p30)$p, shapiro.test(p50)$p,
              shapiro.test(p60)$p, shapiro.test(p100)$p, shapiro.test(p200)$p,
              shapiro.test(p500)$p)
lillie <- c(lillie.test(p5)$p, lillie.test(p10)$p, lillie.test(p15)$p,
              lillie.test(p20)$p, lillie.test(p30)$p, lillie.test(p50)$p,
              lillie.test(p60)$p, lillie.test(p100)$p, lillie.test(p200)$p,
              lillie.test(p500)$p)
pvalue <- data.frame(shapiro, lillie)
rownames(pvalue) <- c("n5", "n10", "n15", "n20", "n30", "n50", "n60",
                        "n100", "n200", "n500")

library(formattable)
formattable(pvalue, list(
  shapiro = formatter("span",
    style = x ~ style(color = ifelse(pvalue$shapiro < 0.05, "red", "green"))),
  lillie = formatter("span",
    style = x ~ style(color = ifelse(pvalue$lillie < 0.05, "red", "green")))))

	shapiro	lillie
n5	1.2e-26	1.3e-144
n10	3.5e-20	4.3e-95
n15	5.9e-16	6.5e-50
n20	1.8e-12	5.0e-33
n30	8.3e-10	1.1e-23
n50	9.0e-07	9.4e-14
n60	2.0e-06	1.2e-15
n100	1.8e-04	1.5e-09
n200	4.8e-04	1.6e-06
n500	2.3e-02	2.9e-04

90% Plantas enfermas

P <- 0.9 # parámetro proporción de plantas enfermas
poblacion <- rbinom(N, 1, P) # Población de plantas

p5 <- Estimar500(5)
p10 <- Estimar500(10)
p15 <- Estimar500(15)
p20 <- Estimar500(20)
p30 <- Estimar500(30)
p50 <- Estimar500(50)
p60 <- Estimar500(60)
p100 <- Estimar500(100)
p200 <- Estimar500(200)
p500 <- Estimar500(500)

par(mfrow = c(2, 5))
boxplot(p5, xlab = "n = 5", las = 2)
abline(h = 0.9, col = "red")
boxplot(p10, xlab = "n = 10", las = 2)
abline(h = 0.9, col = "red")
boxplot(p15, xlab = "n = 15", las = 2)
abline(h = 0.9, col = "red")
boxplot(p20, xlab = "n = 20", las = 2)
abline(h = 0.9, col = "red")
boxplot(p30, xlab = "n = 30", las = 2)
abline(h = 0.9, col = "red")
boxplot(p50, xlab = "n = 50", las = 2)
abline(h = 0.9, col = "red")
boxplot(p60, xlab = "n = 60", las = 2)
abline(h = 0.9, col = "red")
boxplot(p100, xlab = "n = 100", las = 2)
abline(h = 0.9, col = "red")
boxplot(p200, xlab = "n = 200", las = 2)
abline(h = 0.9, col = "red")
boxplot(p500, xlab = "n = 500", las = 2)
abline(h = 0.9, col = "red")

analisisp5 <- rbind(round(mean(p5),3), round(var(p5),3))
analisisp10 <- rbind(round(mean(p10),3), round(var(p10),3))
analisisp15 <- rbind(round(mean(p15),3), round(var(p15),3))
analisisp20 <- rbind(round(mean(p20),3), round(var(p20),3))
analisisp30 <- rbind(round(mean(p30),3), round(var(p30),3))
analisisp50 <- rbind(round(mean(p50),3), round(var(p50),3))
analisisp60 <- rbind(round(mean(p60),3), round(var(p60),3))
analisisp100 <- rbind(round(mean(p100),3), round(var(p100),3))
analisisp200 <- rbind(round(mean(p200),3), round(var(p200),3))
analisisp500 <- rbind(round(mean(p500),3), round(var(p500),3))
analsisisestimadores <- cbind(analisisp5, analisisp10, analisisp15, analisisp20, analisisp30, analisisp50,
analisisp60, analisisp100, analisisp200, analisisp500)
colnames(analsisisestimadores) <- c("n=5", "n=10", "n=15", "n=20", "n=30", "n=50", "n=60", "n=100", "n=200", "n=500")
rownames(analsisisestimadores) <- c("media", "varianza")
print(analsisisestimadores)

           n=5  n=10  n=15  n=20  n=30  n=50  n=60 n=100 n=200 n=500
media    0.907 0.907 0.903 0.903 0.906 0.905 0.902 0.902   0.9   0.9
varianza 0.019 0.008 0.006 0.004 0.003 0.002 0.001 0.001   0.0   0.0

par(mfrow = c(2, 5))
hist(p5)
abline(v = mean(p5), col = "red")
abline(v = median(p5), col = "blue")
hist(p10)
abline(v = mean(p10), col = "red")
abline(v = median(p10), col = "blue")
hist(p15)
abline(v = mean(p15), col = "red")
abline(v = median(p15), col = "blue")
hist(p20)
abline(v = mean(p20), col = "red")
abline(v = median(p20), col = "blue")
hist(p30)
abline(v = mean(p30), col = "red")
abline(v = median(p30), col = "blue")
hist(p50)
abline(v = mean(p50), col = "red")
abline(v = median(p50), col = "blue")
hist(p60)
abline(v = mean(p60), col = "red")
abline(v = median(p60), col = "blue")
hist(p100)
abline(v = mean(p100), col = "red")
abline(v = median(p100), col = "blue")
hist(p200)
abline(v = mean(p200), col = "red")
abline(v = median(p200), col = "blue")
hist(p500)
abline(v = mean(p500), col = "red")
abline(v = median(p500), col = "blue")

par(mfrow = c(2, 5))
qqnorm(p5, main = "Normal Q-Q n=5")
qqline(p5, col = "red")
qqnorm(p10, main = "Normal Q-Q n=10")
qqline(p10, col = "red")
qqnorm(p15, main = "Normal Q-Q n=15")
qqline(p15, col = "red")
qqnorm(p20, main = "Normal Q-Q n=20")
qqline(p20, col = "red")
qqnorm(p30, main = "Normal Q-Q n=30")
qqline(p30, col = "red")
qqnorm(p50, main = "Normal Q-Q  n=50")
qqline(p50, col = "red")
qqnorm(p60, main = "Normal Q-Q n=60")
qqline(p60, col = "red")
qqnorm(p100, main = "Normal Q-Q n=100")
qqline(p100, col = "red")
qqnorm(p200, main = "Normal Q-Q n=200")
qqline(p200, col = "red")
qqnorm(p500, main = "Normal Q-Q n=500")
qqline(p500, col = "red")

library(nortest)
library(formattable)
options(digits = 2)
shapiro <- c(shapiro.test(p5)$p, shapiro.test(p10)$p, shapiro.test(p15)$p,
              shapiro.test(p20)$p, shapiro.test(p30)$p, shapiro.test(p50)$p,
              shapiro.test(p60)$p, shapiro.test(p100)$p, shapiro.test(p200)$p,
              shapiro.test(p500)$p)
lillie <- c(lillie.test(p5)$p, lillie.test(p10)$p, lillie.test(p15)$p,
              lillie.test(p20)$p, lillie.test(p30)$p, lillie.test(p50)$p,
              lillie.test(p60)$p, lillie.test(p100)$p, lillie.test(p200)$p,
              lillie.test(p500)$p)
pvalue <- data.frame(shapiro, lillie)
rownames(pvalue) <- c("n5", "n10", "n15", "n20", "n30", "n50", "n60",
                        "n100", "n200", "n500")

library(formattable)
formattable(pvalue, list(
  shapiro = formatter("span",
    style = x ~ style(color = ifelse(pvalue$shapiro < 0.05, "red", "green"))),
  lillie = formatter("span",
    style = x ~ style(color = ifelse(pvalue$lillie < 0.05, "red", "green")))))

	shapiro	lillie
n5	7.7e-30	1.9e-208
n10	7.4e-23	1.6e-70
n15	2.3e-18	1.8e-63
n20	4.4e-15	2.1e-42
n30	2.2e-12	7.9e-28
n50	1.2e-07	3.6e-19
n60	7.6e-09	1.2e-20
n100	1.7e-05	3.8e-12
n200	3.8e-03	7.4e-06
n500	3.2e-02	1.4e-04

Respecto a las propiedades del estimador, se observa en los 3 casos que no necesariamente se consigue una calidad o efectividad con el aumento en el número de observaciones.

Referente a los histogramas, se observa que a mayor número de muestras tomadas, más simétrica se va volviendo el histograma, lo que corrobora el Teorema del Límite Central, esta misma evidencia, se puede apreciar en los gráficos Q-Q.

Finalmente, las simulaciones con una proporción poblacional de plantas enfermas de 0.1 y 0.9, ninguno de los escenarios simulados fue suficiente para no rechazar la hipóteis de normalidad. Esto se debe seguramente porque, al colocar una proporción poblacional muy baja (0.10), la distribución se sesga a la derecha. De igual manera, cuando la proporción es muy alta (0.9), la distribución se sesga a la izquierda perdiendo simetría.

Problema 4: Estimación Boostrap

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.

Extraer un gran número de muestras (suponga k = 1000). Para cada una de las muestra bootstrap obtenidas calcule la media X∗i

consumo <- c(7.69, 4.97, 4.56, 6.49, 4.34, 6.24, 4.45)
set.seed(123)
muestras <- sample(consumo, 7000, replace = TRUE)
simulacion <- matrix(muestras, nrow = 1000, ncol = 7)
med <- apply(simulacion, 1, mean)

Calcule el intervalo de confianza conformado por los percentiles P2.5 y P97.5

##Método 1: Percentiles P2.5 y P97.5
Met1 <- quantile(med, probs = c(0.025, 0.975))
Met1

 2.5% 97.5% 
  4.7   6.4 

Calcule el intervalo de confianza utilizando el Método 2

#Método 2: Media
x <- mean(med)
lolimit <- as.numeric(2*x - Met1[2])
hilimit <- as.numeric(2*x - Met1[1])
cat("Límite inferior:", lolimit,
"\nLímite superior:", hilimit)

Límite inferior: 4.6 
Límite superior: 6.3

Compare los resultados obtenidos. Comente los resultados. ¿Confiaría en estas estimaciones?

#Grafica de límites
hist(med)
abline(v = Met1, col = "red")
abline(v = lolimit, col = "blue")
abline(v = hilimit, col = "blue")

Se observa que los resultados del intervalo de confianza conseguido por los dos métodos son muy similires con cierto desplazamiento a la izquierda del método 2 (línea azul).

Este método de estimación puede no ser muy preciso, ya que si la muestra no es muy representativa respecto a la población, tal como en el ejercicio anterior, podrían obtenerse resultados simulados en lugar de resultados reales, sin embargo, una de sus ventajas a su vez, es que permite estimar cuando hay muy pocos datos y no se conoce su distribución.

Referencias

Royston, J. P. “Algorithm AS 181: The W Test for Normality.” Journal of the Royal Statistical Society. Series C (Applied Statistics) 31, no. 2 (1982): 176–80. https://doi.org/10.2307/2347986.

Teoría de Probabilidad y Estadística Matemática, Dr. rer. nat. Humberto LLinás Solano, Departamento de Matemáticas y Estadística, Universidad del Norte (Barranquilla, Colombia)