Problema 1

Estimación del valor de \(\pi\)

La siguiente figura sugiere como estimar el valor de π con una simulación. En la figura, un círcuito con un área igual a π/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 é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 la estimación de π/4.De este último resultado se encontrar una aproximación para el valor de π.

Problema tomado de Navidi(2006) 

n <- 10000
x <- runif(n)
y <- runif(n)
d <- (x - 0.5)^2 + (y-0.5)^2
u <- as.numeric(d <= 0.25)
p1 <- sum(u)/n
epi <- (sum(u)/n)*4
error <- pi-epi

cat("proporción de puntos dentro del circulo :", p1, "\n")
cat("estimación de pi                        :", epi, "\n")
cat("ERROR de estimación                     :", error)
proporción de puntos dentro del circulo : 0.7841 
estimación de pi                        : 3.1364 
ERROR de estimación                     : 0.005192654   




library(ggplot2)

n <- 10000

x <- runif(n)
y <- runif(n)

d <- (x - 0.5)^2 + (y - 0.5)^2

dentro <- d <= 0.25

data <- data.frame(x = x, y = y, dentro = dentro)

ggplot(data, aes(x = x, y = y)) +
  geom_point(aes(color = dentro), size = 0.5) +  
  stat_function(fun = function(x) sqrt(0.25 - (x - 0.5)^2) + 0.5, geom = "line", color = "red") +  
  stat_function(fun = function(x) -sqrt(0.25 - (x - 0.5)^2) + 0.5, geom = "line", color = "red") +  
  coord_fixed() +  
  labs(title = "Estimación de Pi", x = "X", y = "Y") +
  theme_minimal()

Se puede concluír que la simulación de monte carlo fue un proceso útil para estimar el valor de pi; esto se evidencia por un error de estimación bajo de 0.005192654 para este caso.

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 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:

θ1ˆ=X1+X26+X3+X43

θ2ˆ=(X1+2X2+3X3+4X4)5

θ3ˆ=X1+X2+X3+X44

θ4ˆ=min{X1,X2,X3,X4}+max{X1,X2,X3,X4}2

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

  2. En cada caso evalue las propiedades de insesgadez, eficiencia y consistencia

  3. Suponga un valor para el parámetro θ

Para el problema 2 se promone un valor Theta = 5. Primero, se genera un análisis con una muestra n= 20

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

data <- data.frame(
theta1 <- (x1 + x2)/6 + (x3 + x4)/3,
theta2 <- (x1 + 2*x2 + 3*x3 + 4*x4)/5,
theta3 <- (x1 + x2 + x3 + x4)/4,
theta4 <- (min(c(x1,x2,x3,x4)) + max(c(x1,x2,x3,x4)))/2
)

names(data) <- c("theta1", "theta2", "theta3", "theta4")

boxplot(data, las=1, col= "#CCFAFF")
abline(h=5, col= "red")

medias <- apply(data, 2, mean)
desv_stds <- apply(data, 2, sd)


tabla_estadisticas <- data.frame(
  Media = medias,
  Desviacion = desv_stds  
)


knitr::kable(tabla_estadisticas, format = "html", caption = "Tabla de Medias y Desviaciones Estándar", table.attr = 'border="1"')
Tabla de Medias y Desviaciones Estándar
Media Desviacion
theta1 5.742275 2.768597
theta2 11.349315 5.878354
theta3 5.691022 2.528215
theta4 10.701018 0.000000



Se evidencia que para para n= 20, Theta sub 1 presenta una media algo aproximada a 5 y una varianza leve lo cual sugiere un estimador insesgado. Theta sub 3 posee una media alejada del valor 5 una alta variabilidad que sugiere un sesgo considerable. Para Theta sub 3 se presenta una media aproximada a 5 y una desviación no tan grande aunque este estimador no presenta resultados tan cercanos como Theta sub 1. Por último Theta 4 posee un sesgo significativo. A continuación se realiza en análisis con una muestra n= 50

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

data <- data.frame(
theta1 <- (x1 + x2)/6 + (x3 + x4)/3,
theta2 <- (x1 + 2*x2 + 3*x3 + 4*x4)/5,
theta3 <- (x1 + x2 + x3 + x4)/4,
theta4 <- (min(c(x1,x2,x3,x4)) + max(c(x1,x2,x3,x4)))/2
)

names(data) <- c("theta1", "theta2", "theta3", "theta4")

boxplot(data, las=1, col= "#CCFAFF")
abline(h=5, col= "red")

medias <- apply(data, 2, mean)
desv_stds <- apply(data, 2, sd)


tabla_estadisticas <- data.frame(
  Media = medias,
  Desviacion = desv_stds  
)


knitr::kable(tabla_estadisticas, format = "html", caption = "Tabla de Medias y Desviaciones Estándar", table.attr = 'border="1"')
Tabla de Medias y Desviaciones Estándar
Media Desviacion
theta1 5.200126 2.568895
theta2 10.491767 5.597645
theta3 5.332894 2.553006
theta4 22.051237 0.000000



Para n= 50 se evidencia que Theta 1 sigue siendo un estimador con una media cercana y una variación baja, para Theta 2 tenemos un sesgo considerable en terminos de la media aunque su varianza no sea tan importante. El estimador Theta 3 ha mejorado y de igual manera sigue siendo uno de los 2 mejores. Theta 4 sigue estando muy por encima del valor objetivo.A continuación se realiza en análisis con una muestra n= 100

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

data <- data.frame(
theta1 <- (x1 + x2)/6 + (x3 + x4)/3,
theta2 <- (x1 + 2*x2 + 3*x3 + 4*x4)/5,
theta3 <- (x1 + x2 + x3 + x4)/4,
theta4 <- (min(c(x1,x2,x3,x4)) + max(c(x1,x2,x3,x4)))/2
)

names(data) <- c("theta1", "theta2", "theta3", "theta4")

boxplot(data, las=1, col= "#CCFAFF")
abline(h=5, col= "red")

medias <- apply(data, 2, mean)
desv_stds <- apply(data, 2, sd)


tabla_estadisticas <- data.frame(
  Media = medias,
  Desviacion = desv_stds  
)


knitr::kable(tabla_estadisticas, format = "html", caption = "Tabla de Medias y Desviaciones Estándar", table.attr = 'border="1"')
Tabla de Medias y Desviaciones Estándar
Media Desviacion
theta1 4.699832 2.440156
theta2 9.609136 5.179683
theta3 4.758946 2.465545
theta4 20.134670 0.000000



Para n= 100 el estimador Theta 2 parece alejarse todavía más del valor objetivo. Los estimadores Theta 1 y Theta 2 parecen alinearse más al valor siendo ahora dos estimadores consistentes. A continuación se realiza en análisis con una muestra n= 1000

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

data <- data.frame(
theta1 <- (x1 + x2)/6 + (x3 + x4)/3,
theta2 <- (x1 + 2*x2 + 3*x3 + 4*x4)/5,
theta3 <- (x1 + x2 + x3 + x4)/4,
theta4 <- (min(c(x1,x2,x3,x4)) + max(c(x1,x2,x3,x4)))/2
)

names(data) <- c("theta1", "theta2", "theta3", "theta4")

boxplot(data, las=1, col= "#CCFAFF")
abline(h=5, col= "red")

medias <- apply(data, 2, mean)
desv_stds <- apply(data, 2, sd)


tabla_estadisticas <- data.frame(
  Media = medias,
  Desviacion = desv_stds  
)


knitr::kable(tabla_estadisticas, format = "html", caption = "Tabla de Medias y Desviaciones Estándar", table.attr = 'border="1"')
Tabla de Medias y Desviaciones Estándar
Media Desviacion
theta1 4.911924 2.577985
theta2 9.833530 5.395270
theta3 4.905984 2.440510
theta4 18.324053 0.000000



Finalmente, para n= 1000 Theta 1 presenta una desviación ligeramente mayor que para n= 100, sin embargo sigue siendo un buen estimador.



Conclusiones generales

Para Theta 1 la consistencia mejora a medida que aumenta el numero n y podemos decir que es insesgado, además, se concluye que es eficiente y consistente debido a que a medida que aumenta el tamaño de la muestra mdisminuye la desviación. Para theta 2 se concluye que no es un estimador con un alto sesgo, para Theta 3 se concluye que es insesgado, eficiente y consistente aunque no tan buen estimador como Theta 1, y para Theta 4 se concluye que es altamente sesgado



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 pˆpara un tamaño de muestra dado n.

Repita el escenario anterior (b) n=500 veces y analice los resultados en cuanto al comportamiento de los 500 resultados del estimador 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.

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.

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.

  set.seed(123)
  n <- 1000
  x <-rbinom(n, 1, 0.5) 

proporcion <- function(m){
  p <- sum(sample(x,m))/m
  return(p)
}
proporcion(10)
par(mfrow = c(1, 3)) 

w <- rep(5,1000)
ensayo1 <- sapply(w, proporcion)

hist(ensayo1, col= "purple", main="Proporciones con muestra 5")

qqnorm(ensayo1, main="Q-Q plot con muestra= 5")  
qqline(ensayo1, col="red")

shapiro.test(ensayo1)
Prueba Shapiro Wilk
  w:0.92989
  p-value < 2.2e-16

par(mfrow = c(1, 2)) 
w <- rep(10,1000)
ensayo2 <- sapply(w, proporcion)

hist(ensayo2, col= "purple", main="Proporciones con muestra 10")

qqnorm(ensayo2, main="Q-Q plot con muestra= 10")  
qqline(ensayo2, col="red")

shapiro.test(ensayo2)
Prueba Shapiro Wilk
  w: 0.92989
  p-value < 2.2e-16

par(mfrow = c(1, 2)) 
w <- rep(15,1000)
ensayo4 <- sapply(w, proporcion)

hist(ensayo4, col= "purple", main="Proporciones con muestra 15")

qqnorm(ensayo4, main="Q-Q plot con muestra= 15")  
qqline(ensayo4, col="red")

shapiro.test(ensayo4)
Prueba Shapiro Wilk
  w: 0.9763
  p-value = 1.057e-11

par(mfrow = c(1, 2)) 
w <- rep(20,1000)
ensayo5 <- sapply(w, proporcion)

hist(ensayo5, col= "purple", main="Proporciones con muestra 20")

qqnorm(ensayo5, main="Q-Q plot con muestra= 20")  
qqline(ensayo5, col="red")

shapiro.test(ensayo5)
Prueba Shapiro Wilk
  w: 0.98315
  p-value = 2.407e-09

par(mfrow = c(1, 2)) 
w <- rep(30,1000)
ensayo6 <- sapply(w, proporcion)

hist(ensayo6, col= "purple", main="Proporciones con muestra 30")

qqnorm(ensayo6, main="Q-Q plot con muestra= 30")  
qqline(ensayo6, col="red")

shapiro.test(ensayo6)
Prueba Shapiro Wilk
  w: 0.98641
  p-value = 5.199e-08

par(mfrow = c(1, 2)) 
w <- rep(50,1000)
ensayo7 <- sapply(w, proporcion)

hist(ensayo7, col= "purple", main="Proporciones con muestra 50")

qqnorm(ensayo7, main="Q-Q plot con muestra= 50")  
qqline(ensayo7, col="red")

shapiro.test(ensayo7)
Prueba Shapiro Wilk
  w: 0.99087
  p-value = 7.234e-06

par(mfrow = c(1, 2)) 
w <- rep(60,1000)
ensayo8 <- sapply(w, proporcion)

hist(ensayo8, col= "purple", main="Proporciones con muestra 60")

qqnorm(ensayo8, main="Q-Q plot con muestra= 60")  
qqline(ensayo8, col="red")

shapiro.test(ensayo8)
Prueba Shapiro Wilk
  w: 0.99283
  p-value = 9.282e-05

par(mfrow = c(1, 2)) 
w <- rep(100,1000)
ensayo9 <- sapply(w, proporcion)

hist(ensayo9, col= "purple", main="Proporciones con muestra 100")

qqnorm(ensayo9, main="Q-Q plot con muestra= 100")  
qqline(ensayo9, col="red")

shapiro.test(ensayo9)
Prueba Shapiro Wilk
  w: 0.99298
  p-value 0.0001144

par(mfrow = c(1, 2)) 
w <- rep(200,1000)
ensayo10 <- sapply(w, proporcion)

hist(ensayo10, col= "purple", main="Proporciones con muestra 200")

qqnorm(ensayo10, main="Q-Q plot con muestra= 200")  
qqline(ensayo10, col="red")

shapiro.test(ensayo10)
Prueba Shapiro Wilk
  w: 0.99696
  p-value = 0.05313

par(mfrow = c(1, 2)) 
w <- rep(500,1000)
ensayo11 <- sapply(w, proporcion)

hist(ensayo11, col= "purple", main="Proporciones con muestra 500")

qqnorm(ensayo11, main="Q-Q plot con muestra= 500")  
qqline(ensayo11, col="red")

shapiro.test(ensayo11)
Prueba Shapiro Wilk
  w: 0.99775
  p-value = 0.1926



Conclusiones: La prueba shapiro wilk al inicio del ensayo con valores pequeños arroja un p-valor inferior a 0.05 lo cual rechaza la afirmación de seguir una distribución normal; sin embargo, a medida que la muestra aumenta se va ajustando cada vez más a una distribución normal. Esto no solo se evidencia con los resultados de la prueba sino que también en las qq-plot se puede observar como cada vez más los datos se van alineando en la diagonal y también en los histogramas se puede observar como los datos se van acercando más unos a otros; haciendo que la distribución se vuelva más simétrica Por último se puede evidenciar que la variabilidad va disminuyendo a medida que la muestra aumenta.



A continuación se realiza el análisis con un 10% de plantas enfermas 

set.seed(123)
n <- 1000
x <- rbinom(n, 1, 0.1)

proporcion <- function(m){
  p <- sum(sample(x,m))/m
  return(p)
}
proporcion(10)
par(mfrow = c(1, 3)) 

w <- rep(5,1000)
ensayo1 <- sapply(w, proporcion)

hist(ensayo1, col= "orange", main="Proporciones con muestra 5")

qqnorm(ensayo1, main="Q-Q plot con muestra 5")  
qqline(ensayo1, col="red")

shapiro.test(ensayo1)
Prueba Shapiro Wilk
  w: 0.92773
  p-value < 2.2e-16

par(mfrow = c(1, 2)) 
w <- rep(10,1000)
ensayo2 <- sapply(w, proporcion)

hist(ensayo2, col= "orange", main="Proporciones con muestra 10")

qqnorm(ensayo2, main="Q-Q plot con muestra 10")  
qqline(ensayo2, col="red")

shapiro.test(ensayo2)
Prueba Shapiro Wilk
  w: 0.96124
  p-value = 1.228e-15

par(mfrow = c(1, 2)) 
w <- rep(15,1000)
ensayo4 <- sapply(w, proporcion)

hist(ensayo4, col= "orange", main="Proporciones con muestra 15")

qqnorm(ensayo4, main="Q-Q plot con muestra 15")  
qqline(ensayo4, col="red")

shapiro.test(ensayo4)
Prueba Shapiro Wilk
  w: 0.97461
  p-value = 3.258e-12

par(mfrow = c(1, 2)) 
w <- rep(20,1000)
ensayo5 <- sapply(w, proporcion)

hist(ensayo5, col= "orange", main="Proporciones con muestra 20")

qqnorm(ensayo5, main="Q-Q plot con muestra 20")  
qqline(ensayo5, col="red")

shapiro.test(ensayo5)
Prueba Shapiro Wilk
  w: 0.97846
  p-value = 5.196e-11

par(mfrow = c(1, 2)) 
w <- rep(30,1000)
ensayo6 <- sapply(w, proporcion)

hist(ensayo6, col= "orange", main="Proporciones con muestra 30")

qqnorm(ensayo6, main="Q-Q plot con muestra 30")  
qqline(ensayo6, col="red")

shapiro.test(ensayo6)
Prueba Shapiro Wilk
  w: 0.98665
  p-value = 6.583e-08

par(mfrow = c(1, 2)) 
w <- rep(50,1000)
ensayo7 <- sapply(w, proporcion)

hist(ensayo7, col= "orange", main="Proporciones con muestra 50")

qqnorm(ensayo7, main="Q-Q plot con muestra 50")  
qqline(ensayo7, col="red")

shapiro.test(ensayo7)
Prueba Shapiro Wilk
  w: 0.98957
  p-value = 1.549e-06

par(mfrow = c(1, 2)) 
w <- rep(60,1000)
ensayo8 <- sapply(w, proporcion)

hist(ensayo8, col= "orange", main="Proporciones con muestra 60")

qqnorm(ensayo8, main="Q-Q plot con muestra 60")  
qqline(ensayo8, col="red")

shapiro.test(ensayo8)
Prueba Shapiro Wilk
  w: 0.99259
  p-value = 6.728e-05

par(mfrow = c(1, 2)) 
w <- rep(100,1000)
ensayo9 <- sapply(w, proporcion)

hist(ensayo9, col= "orange", main="Proporciones con muestra 100")

qqnorm(ensayo9, main="Q-Q plot con muestra 100")  
qqline(ensayo9, col="red")

shapiro.test(ensayo9)
Prueba Shapiro Wilk
  w: 0.99485
  p-value = 0.001778

par(mfrow = c(1, 2)) 
w <- rep(200,1000)
ensayo10 <- sapply(w, proporcion)

hist(ensayo10, col= "orange", main="Proporciones con muestra 200")

qqnorm(ensayo10, main="Q-Q plot con muestra 200")  
qqline(ensayo10, col="red")

shapiro.test(ensayo10)
Prueba Shapiro Wilk
  w: 0.99608
  p-value = 0.01245

par(mfrow = c(1, 2)) 
w <- rep(500,1000)
ensayo11 <- sapply(w, proporcion)

hist(ensayo11, col= "orange", main="Proporciones con muestra 500")

qqnorm(ensayo11, main="Q-Q plot con muestra 500")  
qqline(ensayo11, col="red")

shapiro.test(ensayo11)
Prueba Shapiro Wilk
  w: 0.99715
  p-value = 0.07321



Se evidencia que de con un 10% de plantas enfermas y realizando el mismo procedimiento, también se tiende a aumentar el p-valor, podemos decir entonces que no importa si hay un 10% de plantas enfermas, si el tamaño de la muestra aumenta, los datos tienden a la normalidad. Sin embargo con un n=500 no se llega a superar el p-valor (no llega a ser una distribución normal) esto se debe a que cuando la proporción poblaciónal está más cerca a los extremos la distribución de las proporciones muestrales es más sesgada.



Ahora se realizará el experimento con el 90% de plantas enfermas 

set.seed(123)
n <- 1000
x <- rbinom(n, 1, 0.9)

proporcion <- function(m){
  p <- sum(sample(x,m))/m
  return(p)
}
proporcion(10)
par(mfrow = c(1, 3)) 

w <- rep(5,1000)
ensayo1 <- sapply(w, proporcion)

hist(ensayo1, col= "pink", main="Proporciones con muestra 5")

qqnorm(ensayo1, main="Q-Q plot con muestra 5")  
qqline(ensayo1, col="red")

shapiro.test(ensayo1)
Prueba Shapiro Wilk
  w: 0.92829
  p-value < 2.2e-16

par(mfrow = c(1, 2)) 
w <- rep(10,1000)
ensayo2 <- sapply(w, proporcion)

hist(ensayo2, col= "pink", main="Proporciones con muestra 10")

qqnorm(ensayo2, main="Q-Q plot con muestra 10")  
qqline(ensayo2, col="red")

shapiro.test(ensayo2)
Prueba Shapiro Wilk
  w: 0.96341
  p-value = 3.84e-15

par(mfrow = c(1, 2)) 
w <- rep(15,1000)
ensayo4 <- sapply(w, proporcion)

hist(ensayo4, col= "pink", main="Proporciones con muestra 15")

qqnorm(ensayo4, main="Q-Q plot con muestra 15")  
qqline(ensayo4, col="red")

shapiro.test(ensayo4)
Prueba Shapiro Wilk
  w: 0.976
  p-value = 8.535e-12

par(mfrow = c(1, 2)) 
w <- rep(20,1000)
ensayo5 <- sapply(w, proporcion)

hist(ensayo5, col= "pink", main="Proporciones con muestra 20")

qqnorm(ensayo5, main="Q-Q plot con muestra 20")  
qqline(ensayo5, col="red")

shapiro.test(ensayo5)
Prueba Shapiro Wilk
  w: 0.98128
  p-value = 4.88e-10

par(mfrow = c(1, 2)) 
w <- rep(30,1000)
ensayo6 <- sapply(w, proporcion)

hist(ensayo6, col= "pink", main="Proporciones con muestra 30")

qqnorm(ensayo6, main="Q-Q plot con muestra 30")  
qqline(ensayo6, col="red")

shapiro.test(ensayo6)
## 
##  Shapiro-Wilk normality test
## 
## data:  ensayo6
## W = 0.94792, p-value < 2.2e-16
Prueba Shapiro Wilk
  w: 0.98765
  p-value = 1.846e-07

par(mfrow = c(1, 2)) 
w <- rep(50,1000)
ensayo7 <- sapply(w, proporcion)

hist(ensayo7, col= "pink", main="Proporciones con muestra 50")

qqnorm(ensayo7, main="Q-Q plot con muestra 50")  
qqline(ensayo7, col="red")

shapiro.test(ensayo7)
Prueba Shapiro Wilk
  w: 0.99155
  p-value = 1.713e-05

par(mfrow = c(1, 2)) 
w <- rep(60,1000)
ensayo8 <- sapply(w, proporcion)

hist(ensayo8, col= "pink", main="Proporciones con muestra 60")

qqnorm(ensayo8, main="Q-Q plot con muestra 60")  
qqline(ensayo8, col="red")

shapiro.test(ensayo8)
Prueba Shapiro Wilk
  w: 0.99219
  p-value = 3.934e-05

par(mfrow = c(1, 2)) 
w <- rep(100,1000)
ensayo9 <- sapply(w, proporcion)

hist(ensayo9, col= "pink", main="Proporciones con muestra 100")

qqnorm(ensayo9, main="Q-Q plot con muestra 100")  
qqline(ensayo9, col="red")

shapiro.test(ensayo9)
Prueba Shapiro Wilk
  w: 0.99523
  p-value = 0.003182

par(mfrow = c(1, 2)) 
w <- rep(200,1000)
ensayo10 <- sapply(w, proporcion)

hist(ensayo10, col= "pink", main="Proporciones con muestra 200")

qqnorm(ensayo10, main="Q-Q plot con muestra 200")  
qqline(ensayo10, col="red")

shapiro.test(ensayo10)
Prueba Shapiro Wilk
  w: 0.99632
  p-value = 0.01858

par(mfrow = c(1, 2)) 
w <- rep(500,1000)
ensayo11 <- sapply(w, proporcion)

hist(ensayo11, col= "pink", main="Proporciones con muestra 500")

qqnorm(ensayo11, main="Q-Q plot con muestra 500")  
qqline(ensayo11, col="red")

shapiro.test(ensayo10)
Prueba Shapiro Wilk
  w: 0.99632
  p-value = 0.01858



Se evidencia que de con un 90% de plantas enfermas y realizando el mismo procedimiento, también se tiende a aumentar el p-valor, podemos decir entonces que no importa si hay un 90% de plantas enfermas, si el tamaño de la muestra aumenta, los datos tienden a la normalidad sin embargo al igual que con el ejemplo anterior, con una muestra n= 500 no alcanza la normalidad.

Problema 4

Estimació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 X∗i¯, obteniéndose un valor para cada muestra. El intervalo de confianza queda conformado por los percentiles P2.5 y P97.5. Existen dos métodos para estimarlo:

Método 1 (P2.5;P97.5)

Método 2 (2X¯−P97.5;2X¯−P2.5)

Problema tomado de Navidi(2006)

par(mfrow = c(1, 2)) 
n= 7
replica= 1000

x <-c(7.69, 4.97, 4.56, 6.49, 4.34, 6.24, 4.45)

replica <- sample(x,n*replica, replace = TRUE)

replica <- matrix(replica, ncol = 7)

replicamean <-apply(replica, 1, mean)


q <- quantile(replicamean, c(0.025, 0.975)) 

hist (replicamean, col = "green", main= "Método 1")
abline(v=q, col = "red")


qmetodo2 <- c(2*mean(replicamean)-q[2], 2*mean(replicamean)-q[1])

hist (replicamean, col = "green", main= "Método 2")
abline(v=qmetodo2, col = "red")



A primera vista se puede identificar que no hay una gran diferencia entre los dos métodos; el método 2 es más estrecho que el método 1, se podría concluír que el método 2 es el mejor, sin embargo a continuación se procede a realizar un z-test y un t-test para validar si las diferencias son significativas: 

media_metodo1 <- mean(replicamean)
sd_metodo1 <- sd(replicamean)

replicamean_metodo2 <- 2 * media_metodo1 - replicamean
media_metodo2 <- mean(replicamean_metodo2)
sd_metodo2 <- sd(replicamean_metodo2)

n <- length(replicamean)

z_score <- (media_metodo1 - media_metodo2) / sqrt((sd_metodo1^2 / n) + (sd_metodo2^2 / n))

p_value_z <- 2 * (1 - pnorm(abs(z_score)))

cat("Z-test:\n")
## Z-test:
cat("Z-score:", z_score, "\n")
## Z-score: 0
cat("P-value:", p_value_z, "\n")
## P-value: 1
t_test_result <- t.test(replicamean, replicamean_metodo2, paired = FALSE)

cat("T-test:\n")
## T-test:
print(t_test_result)
## 
##  Welch Two Sample t-test
## 
## data:  replicamean and replicamean_metodo2
## t = 0, df = 1998, p-value = 1
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.03950461  0.03950461
## sample estimates:
## mean of x mean of y 
##  5.541141  5.541141



# Problema 5 ## Relaciones entre la potencia, el tamaño de los efectos y el tamaño de la muestra

El “efecto del tamaño” (o “tamaño del efecto”, en inglés “effect size”) en el contexto de la prueba de hipótesis se refiere a la magnitud de la diferencia o la fuerza de la relación que se está investigando entre las variables. En otras palabras, mide la cantidad de cambio o la importancia práctica de los resultados, más allá de simplemente determinar si una diferencia es estadísticamente significativa. El tamaño del efecto es crucial porque, incluso si una prueba estadística muestra que un resultado es significativo (es decir, rechazas la hipótesis nula), el tamaño del efecto te dice si esa diferencia es realmente importante en un sentido práctico o clínico. Por ejemplo, un estudio podría encontrar que un nuevo medicamento reduce la presión arterial de manera estadísticamente significativa, pero el tamaño del efecto te indicaría si la reducción es lo suficientemente grande como para tener relevancia clínica. En resumen, el tamaño del efecto proporciona una medida complementaria a la significancia estadística, ayudando a interpretar el verdadero impacto o importancia de los resultados encontrados. En este problema, nos centraremos en una aplicación que requiere la aplicación de la prueba t de Student para comparar las medias entre dos grupos. En este contexto evaluaremos cómo el efecto de los tamaños o las diferencias en los tamaños muestrales de los grupos influyen en la potencia de la prueba. De manera formal, la potencia se define como la probabilidad de rechazar la hipótesis nula cuando la hipótesis alternativa es verdadera. De forma más coloquial, la potencia es la capacidad de una prueba estadística para identificar un efecto si este realmente existe. En general, desequilibrios muy marcados en los tamaños de muestra tienden a reducir la potencia estadística, incluso cuando se asocian con tamaños de efecto considerables, lo que aumenta la probabilidad de cometer un error de tipo II. Para fundamentar esta afirmación, debes analizar diferentes resultados computacionales que se presentan a continuación.

Caso 1: Variando los tamaños de los efectos (d) En los códigos del archivo llamado caso1.R, para cada tamaño fijo de los efectos d, se modela la relación entre el tamaño muestral y la potencia (manteniendo constante el nivel de significancia α=0.05). En las figuras se visualizan los resultados para tamaño de efecto muy pequeño (d=0.1), pequeño (d=0.2), mediano (d=0.5) y grande (d=0.8). Repite el análisis usando 5 valores distintos del nivel de significancia. ¿Cambian los resultados? ¿Qué ocurre cuando el tamaño de muestra de los grupos que se comparan es de 20, 60, 100 y 140? Analiza y compara los resultados.

Caso 2: Variando los tamaños muestrales En los códigos del archivo llamado caso2.R, se modela la relación entre el tamaño del efecto d y la potencia (manteniendo constante el nivel de significancia α=0.05). Para ello, se considera los siguientes tamaños de muestra, donde n1 es el número de sujetos en el grupo 1 y n2 es el número de sujetos en el grupo 2:

n1=28, n2=1406: n1representa el 2% del tamaño total de la muestra de 1434.

n1=144, n2=1290: n1 representa el 10 % del tamaño total de la muestra de 1434.

n1=287, n2=1147: n1 representa el 20% del tamaño total de la muestra de 1434.

n1=430, n2=1004: n1representa el 30% del tamaño total de la muestra de 1434.

n1=574, n2=860: n1representa el 40% del tamaño total de la muestra de 1434.

n1=717, n2=717: grupos de igual tamaño (esto es óptimo porque conduce a la potencia más alta para un tamaño de efecto dado).

En la figura resultante, se trazaron las curvas de potencia para la prueba t de Student, en función del tamaño del efecto, asumiendo una tasa de error Tipo I del 5%. La comparación de diferentes curvas de potencia (basadas en el tamaño de la muestra de cada grupo) en el mismo gráfico es una representación visual útil de este análisis. En la figura también se trazó una línea discontinua horizontal en un nivel de potencia aceptable del 80% y una línea vertical en el tamaño del efecto que tendría que estar presente en nuestros datos para alcanzar el 80 % de potencia. Se observa que el tamaño del efecto debe ser superior a 0.54 para alcanzar un nivel de potencia aceptable dados tamaños de grupo altamente desequilibrados de n1=28 y n2=1406, en comparación con todos los demás escenarios que conducen al 100% de potencia. Repite el análisis usando 5 valores distintos del nivel de significancia. ¿Cambian los resultados? ¿Qué ocurre cuando n1=28 y n2=1406? Analiza y compara los resultados.



Referencia: 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) 

d <- seq(0.1, 2, by = 0.1) 
n <- 1:150 


alpha_values <- c(0.01, 0.025, 0.05, 0.1, 0.15)

for (alpha in alpha_values) {
  cat("\n\nNivel de significancia: ", alpha, "\n")
  
  t.test.power.effect <- as.data.frame(do.call("cbind", lapply(1:length(d), function(i) {
    sapply(1:length(n), function(j) {
      power.t.test(n = n[j], d = d[i], sig.level = alpha, power = NULL, type = "two.sample")$power
    })
  })))
  
 
  t.test.power.effect[is.na(t.test.power.effect)] <- 0
  colnames(t.test.power.effect) <- paste(d, "effect size")
  
  
  prueba <- t.test.power.effect 
  

  plot(1, type = "n", frame.plot = FALSE, xlab = "Tamaño muestral", ylab = "Potencia", 
       xlim = c(1, 150), ylim = c(0, 1), main = paste("t-Test con Nivel de Significancia =", alpha), axes = FALSE)
  
  abline(v = seq(0, 150, by = 10), col = "lightgray", lty = "dotted")   
  abline(h = seq(0, 1, by = 0.05), col = "lightgray", lty = "dotted")  
  axis(1, seq(0, 150, by = 10))  
  axis(2, seq(0, 1, by = 0.05)) 
  
  columnas <- 1:ncol(prueba)
  color_linea <- rainbow(length(columnas), alpha = 0.5)
  
  for (i in 1:length(columnas)) {
    lines(1:150, prueba[, columnas[i]], col = color_linea[i], lwd = 3)
  }
}
## 
## 
## Nivel de significancia:  0.01

## 
## 
## Nivel de significancia:  0.025

## 
## 
## Nivel de significancia:  0.05

## 
## 
## Nivel de significancia:  0.1

## 
## 
## Nivel de significancia:  0.15



Análisis: A medida que el tamaño de la muestra aumenta, la potencia también aumenta, para niveles de significancia altos la potencia es más alta