Problema 1

ESTIMACIÓN DE \(\pi\)

La siguiente figura sugiere como estimar el valor de π con una simulación. En la figura, un circuito 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 este, 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 π.

Pasos sugeridos:

A.Genere n coordenadas x: X1, . . . , Xn. 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).

n <- 1000
x = runif(n,0,1)

B.Genere 1000 coordenadas y : Y1,…,Yn, utilizando nuevamente la distribución uniforme con valor mínimo de 0 y valor máximo de 1.

# Definir el número de puntos a generar

n <- 1000

x = runif(n,0,1)

C.Cada punto (Xi,Yi) se encuentra dentro del círculo si su distancia desde el centro (0.5,0.5) es menor a 0.5. Para cada par (Xi,Yi) determine si la distancia desde el centro es menor a 0.5. Esto último se puede realizar al calcular el valor (Xi−0.5)2+(Yi−0.5)2, que es el cuadrado de la distancia, y al determinar si es menor que 0.25.

n <- 1000

# Generar coordenadas aleatorias para x e y
x <- runif(n, 0, 1)
y <- runif(n, 0, 1)

# Calcular la distancia de cada punto al centro (0.5, 0.5)
distancia <- (x - 0.5)^2 + (y - 0.5)^2

# Contar cuántos puntos están dentro del círculo (distancia < radio^2)
contador <- as.numeric(distancia < 0.25)

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

# Contar cuántos puntos están dentro del círculo (distancia < radio^2)
contador <- as.numeric(distancia < 0.25)

cat("Suma de puntos dentro del circulo : ", sum(contador))

Suma de puntos dentro del circulo :  824

# Calcular el área aproximada del círculo usando la proporción de puntos dentro del círculo respecto al total de puntos generados
area_aprox <- sum(contador) / n * 4

#¿Cuál es su estimación de π?



cat("Estimación de pi :", sum(contador)/n*4)

Estimación de pi : 3.296

Conclusión:

Se concluye que entre mas grande el tamaño de la muestra sea , el valor se aproxima a el pi.

n= 100
x = runif(n,0,1)
y = runif(n,0,1)
distancia = (x-0.5)^2+(y-0.5)^2
contador= as.numeric(distancia<0.25)
cat("Numero de puntos ",sum(contador))

Numero de puntos  77

cat("Estimación de pi ",sum(contador)/n*4)

Estimación de pi  3.08

# Definir el número de puntos a generar
n <- 100

# Generar coordenadas aleatorias para x e y
x <- runif(n, 0, 1)
y <- runif(n, 0, 1)

# Calcular la distancia de cada punto al centro (0.5, 0.5)
distancia <- (x - 0.5)^2 + (y - 0.5)^2

# Contar cuántos puntos están dentro del círculo (distancia < radio^2)
contador <- as.numeric(distancia < 0.25)

# Calcular el área aproximada del círculo usando la proporción de puntos dentro del círculo
area_aprox <- sum(contador) / n * 4

# Crear un gráfico de dispersión
plot(x, y, col = ifelse(distancia < 0.25, "#9fc3e7", "#E0E5E9"), pch = 20,
     main = "Simulación de puntos aleatorios en un círculo",
     xlab = "Coordenada X", ylab = "Coordenada Y")

# Dibujar el círculo
theta <- seq(0, 2 * pi, length.out = 100)
circ_x <- 0.5 + 0.5 * cos(theta)
circ_y <- 0.5 + 0.5 * sin(theta)
lines(circ_x, circ_y, col = "#0051a1", lwd = 2)

n= 1000
x = runif(n,0,1)
y = runif(n,0,1)
distancia = (x-0.5)^2+(y-0.5)^2
contador= as.numeric(distancia<0.25)
cat("Numero de puntos ",sum(contador))

Numero de puntos  796

cat("Estimación de pi ",sum(contador)/n*4)

Estimación de pi  3.184

# Definir el número de puntos a generar
n <- 1000

# Generar coordenadas aleatorias para x e y
x <- runif(n, 0, 1)
y <- runif(n, 0, 1)

# Calcular la distancia de cada punto al centro (0.5, 0.5)
distancia <- (x - 0.5)^2 + (y - 0.5)^2

# Contar cuántos puntos están dentro del círculo (distancia < radio^2)
contador <- as.numeric(distancia < 0.25)

# Calcular el área aproximada del círculo usando la proporción de puntos dentro del círculo
area_aprox <- sum(contador) / n * 4

# Crear un gráfico de dispersión
plot(x, y, col = ifelse(distancia < 0.25, "#9fc3e7", "#E0E5E9"), pch = 20,
     main = "Simulación de puntos aleatorios en un círculo",
     xlab = "Coordenada X", ylab = "Coordenada Y")

# Dibujar el círculo
theta <- seq(0, 2 * pi, length.out = 100)
circ_x <- 0.5 + 0.5 * cos(theta)
circ_y <- 0.5 + 0.5 * sin(theta)
lines(circ_x, circ_y, col = "#0051a1", lwd = 2)

n= 10000
x = runif(n,0,1)
y = runif(n,0,1)
distancia = (x-0.5)^2+(y-0.5)^2
contador= as.numeric(distancia<0.25)
cat("Numero de puntos ",sum(contador))

Numero de puntos  7896

cat("Estimación de pi ",sum(contador)/n*4)

Estimación de pi  3.1584

# Definir el número de puntos a generar
n <- 10000

# Generar coordenadas aleatorias para x e y
x <- runif(n, 0, 1)
y <- runif(n, 0, 1)

# Calcular la distancia de cada punto al centro (0.5, 0.5)
distancia <- (x - 0.5)^2 + (y - 0.5)^2

# Contar cuántos puntos están dentro del círculo (distancia < radio^2)
contador <- as.numeric(distancia < 0.25)

# Calcular el área aproximada del círculo usando la proporción de puntos dentro del círculo
area_aprox <- sum(contador) / n * 4

# Crear un gráfico de dispersión
plot(x, y, col = ifelse(distancia < 0.25, "#9fc3e7", "#E0E5E9"), pch = 20,
     main = "Simulación de puntos aleatorios en un círculo",
     xlab = "Coordenada X", ylab = "Coordenada Y")

# Dibujar el círculo
theta <- seq(0, 2 * pi, length.out = 100)
circ_x <- 0.5 + 0.5 * cos(theta)
circ_y <- 0.5 + 0.5 * sin(theta)
lines(circ_x, circ_y, col = "#0051a1", lwd = 2)

n= 100000
x = runif(n,0,1)
y = runif(n,0,1)
distancia = (x-0.5)^2+(y-0.5)^2
contador= as.numeric(distancia<0.25)
cat("Numero de puntos ",sum(contador))

Numero de puntos  78494

cat("Estimación de pi ",sum(contador)/n*4)

Estimación de pi  3.13976

# Definir el número de puntos a generar
n <- 100000

# Generar coordenadas aleatorias para x e y
x <- runif(n, 0, 1)
y <- runif(n, 0, 1)

# Calcular la distancia de cada punto al centro (0.5, 0.5)
distancia <- (x - 0.5)^2 + (y - 0.5)^2

# Contar cuántos puntos están dentro del círculo (distancia < radio^2)
contador <- as.numeric(distancia < 0.25)

# Calcular el área aproximada del círculo usando la proporción de puntos dentro del círculo
area_aprox <- sum(contador) / n * 4

# Crear un gráfico de dispersión
plot(x, y, col = ifelse(distancia < 0.25, "#9fc3e7", "#E0E5E9"), pch = 20,
     main = "Simulación de puntos aleatorios en un círculo",
     xlab = "Coordenada X", ylab = "Coordenada Y")

# Dibujar el círculo
theta <- seq(0, 2 * pi, length.out = 100)
circ_x <- 0.5 + 0.5 * cos(theta)
circ_y <- 0.5 + 0.5 * sin(theta)
lines(circ_x, circ_y, col = "#0051a1", lwd = 2)

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:

\(\hat{\theta}_1 = \frac{X_1 + X_2 + X_3 + X_4}{6}\)
\(\hat{\theta}_2 = \left(\frac{X_1 + 2X_2 + 3X_3 + 4X_4}{5}\right)\)
\(\hat{\theta}_3 = \frac{X_1 + X_2 + X_3 + X_4^4}{4}\)
\(\hat{\theta}_4 = \frac{\min\{X_1, X_2, X_3, X_4\} + \max\{X_1, X_2, X_3, X_4\}}{2}\)

4. Desarrollo

El estimador consiste en un proceso estadístico que intenta predecir el valor de algún parámetro de cierta población, esto con el fin de generalizar el resultado obtenido en la muestra.

Dentro de las propiedades de los estimadores, se encuentran unas características importantes que se deben tener en cuenta para la evaluación, como:

Insesgado: la media de su distribución muestral es igual al parámetro de su población.
Eficiente: el estimador más eficiente es el que tiene la varianza más pequeña.
Consistente: si a medida que n aumenta el valor del estadístico se aproxima al valor del parámetro.
Suficiente: un estimador es suficiente si ningún otro estimador puede proporcionar más información sobre el parámetro.

A continuación se analiza un grupo de estimadores propuestos para la estimación de un parámetro asociado a un modelo de probabilidad, en donde se tienen cuatro muestras siendo n: 20, 50, 100 y 1000. De esta forma, se desea identificar el mejor estimador bajo los parámetros de insesgadez, eficiencia y consistencia.

4.1. Muestra n = 20

# Fijar la semilla aleatoria para reproducibilidad.

set.seed(418) # Se coloca para no tener cambios en el valor en cada ejecución. Se fija la semilla aleatoria.

# Definición de parámetros.

n = 4 # Tamaño de cada muestra.
n_muestras = 20 # Número de cantidades de repetición.
theta <- 12 # Parámetro para la distribución exponencial.

# Generación de datos X a partir de una distribución exponencial.

x = matrix(rexp(n*n_muestras, rate = 1/theta), nrow = n_muestras) # "x = matrix" se crea una matriz llamada "X", con datos generados aleatoriamente de acuerdo a una distribución exponencial. "rexp(n*n_muestras" se genera el número de repeticiones solicitadas de acuerdo al valor otrogado en "n_muestras".

# Definición de los estimadores.

# Se definen cuatro funciones que representan diferentes estimadores.

estimador_1 = function(x) {((x[1] + x [2])/6) + ((x[3] + x [4])/3)}
estimador_2 = function(x) {(x[1] + 2*x[2] + 3*x[3] + 4*x[4] )/5}
estimador_3 = function(x) sum(x)/4 
estimador_4 = function(x) (min(x) + max(x)) / 2

# Cálculo de los estimadores para cada muestra.
# Se crean cuatro matrices llamadas T1, T2, T3, T4; en donde se aplica la función asignada a cada estimador, de acuerdo al punto anterior.
# Con "apply" se indica que se debe aplicar la función a cada matriz generada. "x" es la matriz de la que se desea aplicar la función. "1" indica que se aplicará la función a todas las filas y "estimador" es la función que se aplicará.

T1 = matrix(apply(x, 1, estimador_1), nrow = n_muestras) 
T2 = matrix(apply(x, 1, estimador_2), nrow = n_muestras) 
T3 = matrix(apply(x, 1, estimador_3), nrow = n_muestras) 
T4 = matrix(apply(x, 1, estimador_4), nrow = n_muestras) 
  
# Creación de un dataframe para almacenar los resultados de los estimadores.

T1234 = data.frame(T1, T2, T3, T4) 

# Visualización de los estimadores mediante un diagrama de cajas y bigotes.

boxplot (T1234, las = 1, main = "Comparación estimadores con n = 20")

# Línea que indica el valor de Theta en el diagrama

abline(h = 12, col ="red") #12 equivale al theta digitado al inicio del código.

# Calcular la media de cada estimador

EST_MEDIA_20 <- print(apply (T1234, 2, mean)) # Se realiza variable para llamarlo en el dataframe de la consistencia.

      T1       T2       T3       T4 
11.95154 22.63881 12.94579 16.19289

# Calcular la varianza de cada estimador

EST_VARIANZA_20 <- print(apply (T1234, 2, var)) # Se realiza variable para llamarlo en el dataframe de la eficiencia.

       T1        T2        T3        T4 
 33.95156 121.50936  38.86614  70.53312

Para la muestra 1, siendo n=20, se analiza que el mejor resultado se obtiene con T1, pues su varianza es la más baja, siendo de 33.95 con una media de 11.95; sin embargo resulta ser menos insesgado que T3. Algo similar ocurre con T3, en donde su promedio está muy cerca al valor del parámetro y su varianza también resulta baja con respecto a T2 y T4. Por lo cual se concluye que los grupos con mejor valor estadístico serían T1 y T3.

4.2. Muestra n = 50

# Fijar la semilla aleatoria para reproducibilidad.

set.seed(418) # Se coloca para no tener cambios en el valor en cada ejecución.

# Definición de parámetros.

n = 4 # Tamaño de cada muestra.
n_muestras = 50 # Número de cantidades de repetición.
theta <- 12 # Parámetro para la distribución exponencial.

# Generación de datos X a partir de una distribución exponencial.

x = matrix(rexp(n*n_muestras, rate = 1/theta), nrow = n_muestras)

# Definición de los estimadores.

estimador_1 = function(x) {((x[1] + x [2])/6) + ((x[3] + x [4])/3)}
estimador_2 = function(x) {(x[1] + 2*x[2] + 3*x[3] + 4*x[4] )/5}
estimador_3 = function(x) sum(x)/4 
estimador_4 = function(x) (min(x) + max(x)) / 2

# Cálculo de los estimadores para cada muestra.

T1 = matrix(apply(x, 1, estimador_1), nrow = n_muestras) 
T2 = matrix(apply(x, 1, estimador_2), nrow = n_muestras) 
T3 = matrix(apply(x, 1, estimador_3), nrow = n_muestras) 
T4 = matrix(apply(x, 1, estimador_4), nrow = n_muestras) 
  
# Creación de un dataframe para almacenar los resultados de los estimadores.

T1234 = data.frame(T1, T2, T3, T4) 

# Visualización de los estimadores mediante un diagrama de cajas y bigotes.

boxplot (T1234, las = 1, main = "Comparación estimadores con n = 50")

# Línea que indica el valor de Theta en el diagrama

abline(h = 12, col ="red") #12 equivale al theta digitado al inicio del código.

# Calcular la media de cada estimador

EST_MEDIA_50 <- print(apply (T1234, 2, mean)) # Se realiza variable para llamarlo en el dataframe de la consistencia.

      T1       T2       T3       T4 
13.06600 25.46485 13.05895 15.53341

# Calcular la varianza de cada estimador

EST_VARIANZA_50 <- print(apply (T1234, 2, var)) # Se realiza variable para llamarlo en el dataframe de la eficiencia.

       T1        T2        T3        T4 
 69.87854 330.74026  48.70424 101.84063

Para la muestra 2, en donde n=50, encontramos valores extremos o no esperados en todos los estimadores, sobre todo en T2, en donde su varianza también resulta ser significativamente grande, obteniendo un valor de 300.74. Por otro lado, encontramos nuevamente unos valores más estables en T1 y T3; teniendo una varianza de 69.88 y 48.70 respectivamente; sin embargo, aunque T3 cuenta con un valor extremo menos alejado, es T2 un estimador insesgado, ya que su promedio está demasiado cerca al valor del parámetro.

4.3. Muestra n = 100

# Fijar la semilla aleatoria para reproducibilidad.

set.seed(418) # Se coloca para no tener cambios en el valor en cada ejecución.

# Definición de parámetros.

n = 4 # Tamaño de cada muestra.
n_muestras = 100 # Número de cantidades de repetición.
theta <- 12 # Parámetro para la distribución exponencial.

# Generación de datos X a partir de una distribución exponencial.

x = matrix(rexp(n*n_muestras, rate = 1/theta), nrow = n_muestras)

# Definición de los estimadores.

estimador_1 = function(x) {((x[1] + x [2])/6) + ((x[3] + x [4])/3)}
estimador_2 = function(x) {(x[1] + 2*x[2] + 3*x[3] + 4*x[4] )/5}
estimador_3 = function(x) sum(x)/4 
estimador_4 = function(x) (min(x) + max(x)) / 2

# Cálculo de los estimadores para cada muestra.

T1 = matrix(apply(x, 1, estimador_1), nrow = n_muestras) 
T2 = matrix(apply(x, 1, estimador_2), nrow = n_muestras) 
T3 = matrix(apply(x, 1, estimador_3), nrow = n_muestras) 
T4 = matrix(apply(x, 1, estimador_4), nrow = n_muestras) 
  
# Creación de un dataframe para almacenar los resultados de los estimadores.

T1234 = data.frame(T1, T2, T3, T4) 

# Visualización de los estimadores mediante un diagrama de cajas y bigotes.

boxplot (T1234, las = 1, main = "Comparación estimadores con n = 100")

# Línea que indica el valor de Theta en el diagrama

abline(h = 12, col ="red") #12 equivale al theta digitado al inicio del código.

# Calcular la media de cada estimador

EST_MEDIA_100 <- print(apply (T1234, 2, mean)) # Se realiza variable para llamarlo en el dataframe de la consistencia.

      T1       T2       T3       T4 
12.42388 24.49794 12.58265 14.84110

# Calcular la varianza de cada estimador

EST_VARIANZA_100 <- print(apply (T1234, 2, var)) # Se realiza variable para llamarlo en el dataframe de la eficiencia.

       T1        T2        T3        T4 
 41.28078 172.11229  40.42906  84.18982

Para la muestra 3, Cuando n=100, los estimadores T1, T3 y T4 resultan ser más consistentes entre ellos, dejando a T2 nuevamente con la varianza más grande y con una mediana muy por encima a la media total de la población; por otro lado, su mediana es significativamente alta en comparación con el resto de los estimadores.

T1 y T3 nuevamente resultan ser los mejores estimadores, con una varianza de 41.28 y 40.43 respectivamente; además sus medianas se encuentran muy cerca a la media, aunque para el caso de T1 cuenta con valores extremos más alejados que T3.

4.4. Muestra n = 1000

# Fijar la semilla aleatoria para reproducibilidad.

set.seed(418) # Se coloca para no tener cambios en el valor en cada ejecución.

# Definición de parámetros.

n = 4 # Tamaño de cada muestra.
n_muestras = 1000 # Número de cantidades de repetición.
theta <- 12 # Parámetro para la distribución exponencial.

# Generación de datos X a partir de una distribución exponencial.

x = matrix(rexp(n*n_muestras, rate = 1/theta), nrow = n_muestras)

# Definición de los estimadores.

estimador_1 = function(x) {((x[1] + x [2])/6) + ((x[3] + x [4])/3)}
estimador_2 = function(x) {(x[1] + 2*x[2] + 3*x[3] + 4*x[4] )/5}
estimador_3 = function(x) sum(x)/4 
estimador_4 = function(x) (min(x) + max(x)) / 2

# Cálculo de los estimadores para cada muestra.

T1 = matrix(apply(x, 1, estimador_1), nrow = n_muestras) 
T2 = matrix(apply(x, 1, estimador_2), nrow = n_muestras) 
T3 = matrix(apply(x, 1, estimador_3), nrow = n_muestras) 
T4 = matrix(apply(x, 1, estimador_4), nrow = n_muestras) 
  
# Creación de un dataframe para almacenar los resultados de los estimadores.

T1234 = data.frame(T1, T2, T3, T4) 

# Visualización de los estimadores mediante un diagrama de cajas y bigotes.

boxplot (T1234, las = 1, main = "Comparación estimadores con n = 1000")

# Línea que indica el valor de Theta en el diagrama

abline(h = 12, col ="red") #12 equivale al theta digitado al inicio del código.

# Calcular la media de cada estimador

EST_MEDIA_1000 <- print(apply (T1234, 2, mean)) # Se realiza variable para llamarlo en el dataframe de la consistencia.

      T1       T2       T3       T4 
12.00459 24.03566 12.01127 13.89701

# Calcular la varianza de cada estimador

EST_VARIANZA_1000 <- print(apply (T1234, 2, var)) # Se realiza variable para llamarlo en el dataframe de la eficiencia.

       T1        T2        T3        T4 
 43.80968 190.01993  39.49763  58.18627

Cuando n=1000 en la muestra 4, se evidencian valores extremos en todos los estimadores. T1, T3 y T4 resultan tener los mejores resultados. Aunque T4 tiene la varianza más grande de los otros dos estimadores en comparación, siendo esta de 58.18, se encuentra mucho más cerca a la media con valores más extremos.

4. Consistencia

# Encontrar la media

media_identificada <- data.frame(
  muestra_20 = tail(EST_MEDIA_20, 4),
  muestra_50 = tail(EST_MEDIA_50, 4),
  muestra_100 = tail(EST_MEDIA_100, 4),
  muestra_1000 = tail(EST_MEDIA_1000, 4)
)

media_identificada

   muestra_20 muestra_50 muestra_100 muestra_1000
T1   11.95154   13.06600    12.42388     12.00459
T2   22.63881   25.46485    24.49794     24.03566
T3   12.94579   13.05895    12.58265     12.01127
T4   16.19289   15.53341    14.84110     13.89701

La muestra 4 cuando n=1000, se visualiza una mejor consistencia para los estimadores 1 y 3. El estimador 2 no genera un resultado consistente.

4. Sesgo

preal <- 12

tamaños <- c(20, 50, 100, 1000)
sesgo_resultados <- data.frame(n = integer(), estimador = character(), sesgo = numeric())

for(i in tamaños) {
  x = matrix(rexp(n*i, rate = 1/theta), nrow = i)
  T1 = apply(x, 1, estimador_1)
  T2 = apply(x, 1, estimador_2) 
  T3 = apply(x, 1, estimador_3) 
  T4 = apply(x, 1, estimador_4)
  
  sesgo_T1 <- mean(T1) - preal
  sesgo_T2 <- mean(T2) - preal
  sesgo_T3 <- mean(T3) - preal
  sesgo_T4 <- mean(T4) - preal

  sesgo_resultados <- rbind(sesgo_resultados, data.frame(
    n=rep(i, 4), estimador = c("T1", "T2", "T3", "T4"), sesgo = c(sesgo_T1, sesgo_T2, sesgo_T3, sesgo_T4)))
  
  
}

print(sesgo_resultados)

      n estimador      sesgo
1    20        T1 -1.8630568
2    20        T2  7.5198524
3    20        T3 -1.6119687
4    20        T4 -0.2965256
5    50        T1 -0.8044428
6    50        T2 10.2029037
7    50        T3 -0.6101756
8    50        T4  1.3178604
9   100        T1  1.1258249
10  100        T2 14.1890950
11  100        T3  1.0163688
12  100        T4  2.6993690
13 1000        T1  0.4044995
14 1000        T2 12.9164385
15 1000        T3  0.3035003
16 1000        T4  2.2384681

El resultado más insesgado resulta ser en el estimador 3 cuando n=1000. Es importante mencionar que el estimador 1 también cuenta con un mejor sesgo cuando la muestra incrementa. Para efectos del estimador 2, resulta ser sesgado.

4. Eficiencia

# Encontrar varianza más baja de los estimadores

variacion <- data.frame(
  muestra_20 = tail(EST_VARIANZA_20, 4),
  muestra_50 = tail(EST_VARIANZA_50, 4),
  muestra_100 = tail(EST_VARIANZA_100, 4),
  muestra_1000 = tail(EST_VARIANZA_1000, 4)
)

variacion

   muestra_20 muestra_50 muestra_100 muestra_1000
T1   33.95156   69.87854    41.28078     43.80968
T2  121.50936  330.74026   172.11229    190.01993
T3   38.86614   48.70424    40.42906     39.49763
T4   70.53312  101.84063    84.18982     58.18627

La mejor varianza se obtiene cuando es n = 1000 en el estimador 3.

4.5. Conclusión - Problema 2

Cuando se tiene mayor repeticiones dentro de una muestra, se evidencia una mejor insesgadez, eficiencia y consistencia para todos los estimadores; sin embargo, T2 debe ser retirado del análisis, ya que representa mayor desviación y se encuentra más alejado a la media.

Con el incremento de las muestras se evidencia que T4 está cada vez más cerca de la media y con un mejor resultado en su varianza.

Cuando n=1000 se encuentran datos consistentes para T1, T3 y T4.

Se evidencia un problema en todos los estimadores, puesto que los valores extremos o no esperados superan las cifras matemáticamente permitidos para estimaciones, lo que sería necesario realizar una investigación para determinar si hubo un error durante el proceso de toma de la muestra.

El estimador 3, cuando n=1000, resulta ser el más óptimo.

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:

A.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%.

poblacion = c(rep(1,500),rep(0,500))
barplot(table(poblacion),las=1, col= c("#e0e5e9","#9fc3e7"))

B.Genere una función que permita:

Obtener una muestra aleatoria de la población y

#Genere una función que permita obtener una muestra aleatoria de la población y 


muestra = function(x,n){sample(x,n,replace = TRUE)}

muestra = function(x,n){sample(x,n,replace = TRUE)}

Calcule el estimador de la proporción muestral \(\widehat{p}\) para un tamaño de muestra dado n.
- tamaño de las muestras : n = 5
- número de muestras : m = 10

n=5
m=10

y= matrix(muestra(poblacion,n*m), ncol=n)

#y

Calcule el estimador de la proporción muestral para un tamaño de muestra dado n.

#calcule el estimador de la proporción muestral para un tamaño de muestra dado n.

p=apply(y,1,sum)/n

hist(p, las=1, main = "proporción muestral", col = "#9fc3e7")

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

poblacion = c(rep(1,500),rep(0,500))

n = 500

muestra = function(x,n){sample(x,n,replace = TRUE)}
n=500
m=500

y= matrix(muestra(poblacion,n*m), ncol=m)

#y

p=apply(y,1,sum)/n

hist(p,main = "Proporción, n=500",ylab = "Frecuencia", xlab = "p", col="#9fc3e7")

summarytools::descr(p)

Descriptive Statistics  
p  
N: 500  

                         p
----------------- --------
             Mean     0.50
          Std.Dev     0.02
              Min     0.44
               Q1     0.48
           Median     0.50
               Q3     0.52
              Max     0.56
              MAD     0.02
              IQR     0.03
               CV     0.04
         Skewness    -0.12
      SE.Skewness     0.11
         Kurtosis    -0.42
          N.Valid   500.00
        Pct.Valid   100.00

Dado que el coeficiente de asimetria (skewness) es igual a 0.08 la distribución de p un tamaño de muestra n=500 es simetrico( simetrico es cuando el valor es cercano a cero o igual). En cuanto a la variabilidad se puede medir con el coeficiente de variación (CV) que tiene un valor de 0.05, lo que indica que los resultados son muy homogeneos.

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

# n=5

poblacion = c(rep(1,500),rep(0,500))
muestra = function(x,n){sample(x,n,replace = TRUE)}
n=5
m=500

y5= matrix(muestra(poblacion,n*m), ncol=n)


p5=apply(y5,1,sum)/n

par(mfrow = c(1, 2))
hist(p5,main = "Proporción, n=5",ylab = "Frecuencia", xlab = "p", col="#9fc3e7")
qqnorm(p5,main = " ", ylab = "Percenties Muestra", xlab = "Percentiles Normal", col="#0051a1")

S5=shapiro.test(p5)
S5


    Shapiro-Wilk normality test

data:  p5
W = 0.93024, p-value = 1.656e-14

#n=10

poblacion = c(rep(1,500),rep(0,500))
muestra = function(x,n){sample(x,n,replace = TRUE)}
n=10
m=500

y5= matrix(muestra(poblacion,n*m), ncol=n)


p6=apply(y5,1,sum)/n

par(mfrow = c(1, 2))
hist(p6,main = "Proporción, n=10",ylab = "Frecuencia", xlab = "p", col="#9fc3e7")
qqnorm(p6,main = " ", ylab = "Percenties Muestra", xlab = "Percentiles Normal", col="#0051a1")

S6=shapiro.test(p6)
S6


    Shapiro-Wilk normality test

data:  p6
W = 0.96738, p-value = 4.24e-09

#n=15

poblacion = c(rep(1,500),rep(0,500))
muestra = function(x,n){sample(x,n,replace = TRUE)}
n=15
m=500

y5= matrix(muestra(poblacion,n*m), ncol=n)


p7=apply(y5,1,sum)/n

par(mfrow = c(1, 2))
hist(p7,main = "Proporción, n=15",ylab = "Frecuencia", xlab = "p", col="#9fc3e7")
qqnorm(p7,main = " ", ylab = "Percenties Muestra", xlab = "Percentiles Normal", col="#0051a1")

S7=shapiro.test(p7)
S7


    Shapiro-Wilk normality test

data:  p7
W = 0.97244, p-value = 4.328e-08

#n=20

poblacion = c(rep(1,500),rep(0,500))
muestra = function(x,n){sample(x,n,replace = TRUE)}
n=20
m=500

y5= matrix(muestra(poblacion,n*m), ncol=n)


p8=apply(y5,1,sum)/n

par(mfrow = c(1, 2))
hist(p8,main = "Proporción, n=20",ylab = "Frecuencia", xlab = "p", col="#9fc3e7")
qqnorm(p8,main = " ", ylab = "Percenties Muestra", xlab = "Percentiles Normal", col="#0051a1")

S8=shapiro.test(p8)
S8


    Shapiro-Wilk normality test

data:  p8
W = 0.98046, p-value = 3.004e-06

#n=30


poblacion = c(rep(1,500),rep(0,500))
muestra = function(x,n){sample(x,n,replace = TRUE)}
n=30
m=500

y5= matrix(muestra(poblacion,n*m), ncol=n)


p9=apply(y5,1,sum)/n

par(mfrow = c(1, 2))
hist(p9,main = "Proporción, n=30",ylab = "Frecuencia", xlab = "p", col="#9fc3e7")
qqnorm(p9,main = " ", ylab = "Percenties Muestra", xlab = "Percentiles Normal", col="#0051a1")

S9=shapiro.test(p9)
S9


    Shapiro-Wilk normality test

data:  p9
W = 0.98638, p-value = 0.0001266

#n=50

poblacion = c(rep(1,500),rep(0,500))
muestra = function(x,n){sample(x,n,replace = TRUE)}
n=50
m=500

y5= matrix(muestra(poblacion,n*m), ncol=n)


p10=apply(y5,1,sum)/n

par(mfrow = c(1, 2))
hist(p10,main = "Proporción, n=50",ylab = "Frecuencia", xlab = "p", col="#9fc3e7")
qqnorm(p10,main = " ", ylab = "Percenties Muestra", xlab = "Percentiles Normal", col="#0051a1")

S10=shapiro.test(p10)
S10


    Shapiro-Wilk normality test

data:  p10
W = 0.98783, p-value = 0.0003519

#n=60

poblacion = c(rep(1,500),rep(0,500))
muestra = function(x,n){sample(x,n,replace = TRUE)}
n=60
m=500

y5= matrix(muestra(poblacion,n*m), ncol=n)


p11=apply(y5,1,sum)/n

par(mfrow = c(1, 2))
hist(p11,main = "Proporción, n=60",ylab = "Frecuencia", xlab = "p", col="#9fc3e7")
qqnorm(p11,main = " ", ylab = "Percenties Muestra", xlab = "Percentiles Normal", col="#0051a1")

S11=shapiro.test(p11)
S11


    Shapiro-Wilk normality test

data:  p11
W = 0.99307, p-value = 0.02097

#n=100

poblacion = c(rep(1,500),rep(0,500))
muestra = function(x,n){sample(x,n,replace = TRUE)}
n=100
m=500

y5= matrix(muestra(poblacion,n*m), ncol=n)


p12=apply(y5,1,sum)/n

par(mfrow = c(1, 2))
hist(p12,main = "Proporción, n=100",ylab = "Frecuencia", xlab = "p", col="#9fc3e7")
qqnorm(p12,main = " ", ylab = "Percenties Muestra", xlab = "Percentiles Normal", col="#0051a1")

S12=shapiro.test(p12)
S12


    Shapiro-Wilk normality test

data:  p12
W = 0.99377, p-value = 0.03755

#n=200

poblacion = c(rep(1,500),rep(0,500))
muestra = function(x,n){sample(x,n,replace = TRUE)}
n=200
m=500

y5= matrix(muestra(poblacion,n*m), ncol=n)


p13=apply(y5,1,sum)/n

par(mfrow = c(1, 2))
hist(p13,main = "Proporción, n=200",ylab = "Frecuencia", xlab = "p", col="#9fc3e7")
qqnorm(p13,main = " ", ylab = "Percenties Muestra", xlab = "Percentiles Normal", col="#0051a1")

S13=shapiro.test(p13)
S13


    Shapiro-Wilk normality test

data:  p13
W = 0.99608, p-value = 0.2521

#n=500

poblacion = c(rep(1,500),rep(0,500))
muestra = function(x,n){sample(x,n,replace = TRUE)}
n=500
m=500

y5= matrix(muestra(poblacion,n*m), ncol=n)


p14=apply(y5,1,sum)/n

par(mfrow = c(1, 2))
hist(p14,main = "Proporción, n=500",ylab = "Frecuencia", xlab = "p", col="#9fc3e7")
qqnorm(p14,main = " ", ylab = "Percenties Muestra", xlab = "Percentiles Normal", col="#0051a1")

S14=shapiro.test(p14)
S14


    Shapiro-Wilk normality test

data:  p14
W = 0.99626, p-value = 0.291

library(knitr)
library(kableExtra)
library(tidyverse)
# Crear la tabla
d <- data.frame(
  n = c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500),
  valor.p = c(S5$p.value, S6$p.value, S7$p.value, S8$p.value, S9$p.value, S10$p.value, S11$p.value, S12$p.value, S13$p.value, S14$p.value),
  normalidad = c("no", "no", "no", "no", "no", "no", "no", "no", "si", "si")
)

# Crear la tabla con kable y kableExtra
styled_table <- d %>%
  mutate(
    normalidad = ifelse(normalidad == "si", "Sí", "No")
  ) %>%
  kable("html") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE) %>%
  column_spec(1, width = "70px") %>%
  column_spec(2, width = "100px") %>%
  column_spec(3, width = "100px")

styled_table

n	valor.p	normalidad
5	0.0000000	No
10	0.0000000	No
15	0.0000000	No
20	0.0000030	No
30	0.0001266	No
50	0.0003519	No
60	0.0209676	No
100	0.0375501	No
200	0.2520881	Sí
500	0.2909714	Sí

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

A.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 10%.

poblacion = c(rep(1,100),rep(0,900))
barplot(table(poblacion),las=1, col= c("#e0e5e9","#9fc3e7"))

B.Genere una función que permita:

Obtener una muestra aleatoria de la población y

#Genere una función que permita obtener una muestra aleatoria de la población y 


muestra = function(x,n){sample(x,n,replace = TRUE)}

muestra = function(x,n){sample(x,n,replace = TRUE)}

Calcule el estimador de la proporción muestral \(\widehat{p}\) para un tamaño de muestra dado n.
- tamaño de las muestras : n = 5
- número de muestras : m = 10

n=5
m=10

y= matrix(muestra(poblacion,n*m), ncol=n)

#y

Calcule el estimador de la proporción muestral para un tamaño de muestra dado n.

#calcule el estimador de la proporción muestral para un tamaño de muestra dado n.

p=apply(y,1,sum)/n

hist(p, las=1, main = "proporción muestral", col = "#9fc3e7")

poblacion = c(rep(1,100),rep(0,900))

n = 500

muestra = function(x,n){sample(x,n,replace = TRUE)}
n=500
m=500

y= matrix(muestra(poblacion,n*m), ncol=m)

#y

p=apply(y,1,sum)/n

hist(p,main = "Proporción, n=500",ylab = "Frecuencia", xlab = "p", col="#9fc3e7")

# n=5

poblacion = c(rep(1,100),rep(0,900))
muestra = function(x,n){sample(x,n,replace = TRUE)}
n=5
m=500

y5= matrix(muestra(poblacion,n*m), ncol=n)


p5=apply(y5,1,sum)/n

par(mfrow = c(1, 2))
hist(p5,main = "Proporción, n=5",ylab = "Frecuencia", xlab = "p", col="#9fc3e7")
qqnorm(p5,main = " ", ylab = "Percenties Muestra", xlab = "Percentiles Normal", col="#0051a1")

S5=shapiro.test(p5)
S5


    Shapiro-Wilk normality test

data:  p5
W = 0.72848, p-value < 2.2e-16

#n=10

poblacion = c(rep(1,100),rep(0,900))
muestra = function(x,n){sample(x,n,replace = TRUE)}
n=10
m=500

y5= matrix(muestra(poblacion,n*m), ncol=n)


p6=apply(y5,1,sum)/n

par(mfrow = c(1, 2))
hist(p6,main = "Proporción, n=10",ylab = "Frecuencia", xlab = "p", col="#9fc3e7")
qqnorm(p6,main = " ", ylab = "Percenties Muestra", xlab = "Percentiles Normal", col="#0051a1")

S6=shapiro.test(p6)
S6


    Shapiro-Wilk normality test

data:  p6
W = 0.83392, p-value < 2.2e-16

#n=15

poblacion = c(rep(1,100),rep(0,900))
muestra = function(x,n){sample(x,n,replace = TRUE)}
n=15
m=500

y5= matrix(muestra(poblacion,n*m), ncol=n)


p7=apply(y5,1,sum)/n

par(mfrow = c(1, 2))
hist(p7,main = "Proporción, n=15",ylab = "Frecuencia", xlab = "p", col="#9fc3e7")
qqnorm(p7,main = " ", ylab = "Percenties Muestra", xlab = "Percentiles Normal", col="#0051a1")

S7=shapiro.test(p7)
S7


    Shapiro-Wilk normality test

data:  p7
W = 0.9055, p-value < 2.2e-16

#n=20

poblacion = c(rep(1,100),rep(0,900))
muestra = function(x,n){sample(x,n,replace = TRUE)}
n=20
m=500

y5= matrix(muestra(poblacion,n*m), ncol=n)


p8=apply(y5,1,sum)/n

par(mfrow = c(1, 2))
hist(p8,main = "Proporción, n=20",ylab = "Frecuencia", xlab = "p", col="#9fc3e7")
qqnorm(p8,main = " ", ylab = "Percenties Muestra", xlab = "Percentiles Normal", col="#0051a1")

S8=shapiro.test(p8)
S8


    Shapiro-Wilk normality test

data:  p8
W = 0.92528, p-value = 4.571e-15

#n=30


poblacion = c(rep(1,100),rep(0,900))
muestra = function(x,n){sample(x,n,replace = TRUE)}
n=30
m=500

y5= matrix(muestra(poblacion,n*m), ncol=n)


p9=apply(y5,1,sum)/n

par(mfrow = c(1, 2))
hist(p9,main = "Proporción, n=30",ylab = "Frecuencia", xlab = "p", col="#9fc3e7")
qqnorm(p9,main = " ", ylab = "Percenties Muestra", xlab = "Percentiles Normal", col="#0051a1")

S9=shapiro.test(p9)
S9


    Shapiro-Wilk normality test

data:  p9
W = 0.95672, p-value = 6.054e-11

#n=50

poblacion = c(rep(1,100),rep(0,900))
muestra = function(x,n){sample(x,n,replace = TRUE)}
n=50
m=500

y5= matrix(muestra(poblacion,n*m), ncol=n)


p10=apply(y5,1,sum)/n

par(mfrow = c(1, 2))
hist(p10,main = "Proporción, n=50",ylab = "Frecuencia", xlab = "p", col="#9fc3e7")
qqnorm(p10,main = " ", ylab = "Percenties Muestra", xlab = "Percentiles Normal", col="#0051a1")

S10=shapiro.test(p10)
S10


    Shapiro-Wilk normality test

data:  p10
W = 0.97081, p-value = 1.994e-08

#n=60

poblacion = c(rep(1,100),rep(0,900))
muestra = function(x,n){sample(x,n,replace = TRUE)}
n=60
m=500

y5= matrix(muestra(poblacion,n*m), ncol=n)


p11=apply(y5,1,sum)/n

par(mfrow = c(1, 2))
hist(p11,main = "Proporción, n=60",ylab = "Frecuencia", xlab = "p", col="#9fc3e7")
qqnorm(p11,main = " ", ylab = "Percenties Muestra", xlab = "Percentiles Normal", col="#0051a1")

S11=shapiro.test(p11)
S11


    Shapiro-Wilk normality test

data:  p11
W = 0.96651, p-value = 2.92e-09

#n=100

poblacion = c(rep(1,100),rep(0,900))
muestra = function(x,n){sample(x,n,replace = TRUE)}
n=100
m=500

y5= matrix(muestra(poblacion,n*m), ncol=n)


p12=apply(y5,1,sum)/n

par(mfrow = c(1, 2))
hist(p12,main = "Proporción, n=100",ylab = "Frecuencia", xlab = "p", col="#9fc3e7")
qqnorm(p12,main = " ", ylab = "Percenties Muestra", xlab = "Percentiles Normal", col="#0051a1")

S12=shapiro.test(p12)
S12


    Shapiro-Wilk normality test

data:  p12
W = 0.97913, p-value = 1.402e-06

#n=200

poblacion = c(rep(1,100),rep(0,900))
muestra = function(x,n){sample(x,n,replace = TRUE)}
n=200
m=500

y5= matrix(muestra(poblacion,n*m), ncol=n)


p13=apply(y5,1,sum)/n

par(mfrow = c(1, 2))
hist(p13,main = "Proporción, n=200",ylab = "Frecuencia", xlab = "p", col="#9fc3e7")
qqnorm(p13,main = " ", ylab = "Percenties Muestra", xlab = "Percentiles Normal", col="#0051a1")

S13=shapiro.test(p13)
S13


    Shapiro-Wilk normality test

data:  p13
W = 0.98644, p-value = 0.0001322

#n=500

poblacion = c(rep(1,100),rep(0,900))
muestra = function(x,n){sample(x,n,replace = TRUE)}
n=500
m=500

y5= matrix(muestra(poblacion,n*m), ncol=n)


p14=apply(y5,1,sum)/n

par(mfrow = c(1, 2))
hist(p14,main = "Proporción, n=500",ylab = "Frecuencia", xlab = "p", col="#9fc3e7")
qqnorm(p14,main = " ", ylab = "Percenties Muestra", xlab = "Percentiles Normal", col="#0051a1")

S14=shapiro.test(p14)
S14


    Shapiro-Wilk normality test

data:  p14
W = 0.99153, p-value = 0.005926

#n=1000

poblacion = c(rep(1,100),rep(0,900))
muestra = function(x,n){sample(x,n,replace = TRUE)}
n=1000
m=500

y5= matrix(muestra(poblacion,n*m), ncol=n)


p15=apply(y5,1,sum)/n

par(mfrow = c(1, 2))
hist(p15,main = "Proporción, n=1000",ylab = "Frecuencia", xlab = "p", col="#9fc3e7")
qqnorm(p15,main = " ", ylab = "Percenties Muestra", xlab = "Percentiles Normal", col="#0051a1")

S15=shapiro.test(p15)
S15


    Shapiro-Wilk normality test

data:  p15
W = 0.99709, p-value = 0.5189

library(knitr)
library(kableExtra)
library(tidyverse)
# Crear la tabla
d <- data.frame(
  n = c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500, 1000),
  valor.p = c(S5$p.value, S6$p.value, S7$p.value, S8$p.value, S9$p.value, S10$p.value, S11$p.value, S12$p.value, S13$p.value, S14$p.value, S15$p.value),
  normalidad = c("no", "no", "no", "no", "no", "no", "no", "no", "no", "no","si")
)

# Crear la tabla con kable y kableExtra
styled_table <- d %>%
  mutate(
    normalidad = ifelse(normalidad == "si", "Sí", "No")
  ) %>%
  kable("html") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE) %>%
  column_spec(1, width = "70px") %>%
  column_spec(2, width = "100px") %>%
  column_spec(3, width = "100px")

styled_table

n	valor.p	normalidad
5	0.0000000	No
10	0.0000000	No
15	0.0000000	No
20	0.0000000	No
30	0.0000000	No
50	0.0000000	No
60	0.0000000	No
100	0.0000014	No
200	0.0001322	No
500	0.0059261	No
1000	0.5188586	Sí

En las graficas como en la tabla anterior se puede concluir que en la prueba de shapiro cuando el valor de p es mayor a 0.05 no se rechaza la hipotesis de normalidad, tambien en la tabla muestra que las proporciones para todos los n ninguno es normal, pero decidi agregarun n=1000 que si es normal, para esto se necesitan un mayor tamaño de muestra para encontrar la normalidad.

A.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 90%.

poblacion = c(rep(1,900),rep(0,100))
barplot(table(poblacion),las=1, col= c("#e0e5e9","#9fc3e7"))

B.Genere una función que permita:

Obtener una muestra aleatoria de la población y

#Genere una función que permita obtener una muestra aleatoria de la población y 


muestra = function(x,n){sample(x,n,replace = TRUE)}

muestra = function(x,n){sample(x,n,replace = TRUE)}

Calcule el estimador de la proporción muestral \(\widehat{p}\) para un tamaño de muestra dado n.
- tamaño de las muestras : n = 5
- número de muestras : m = 10

n=5
m=10

y= matrix(muestra(poblacion,n*m), ncol=n)

#y

Calcule el estimador de la proporción muestral para un tamaño de muestra dado n.

#calcule el estimador de la proporción muestral para un tamaño de muestra dado n.

p=apply(y,1,sum)/n

hist(p, las=1, main = "proporción muestral", col = "#9fc3e7")

poblacion = c(rep(1,900),rep(0,100))

n = 500

muestra = function(x,n){sample(x,n,replace = TRUE)}
n=500
m=500

y= matrix(muestra(poblacion,n*m), ncol=m)

#y

p=apply(y,1,sum)/n

hist(p,main = "Proporción, n=500",ylab = "Frecuencia", xlab = "p", col="#9fc3e7")

# n=5

poblacion = c(rep(1,900),rep(0,100))
muestra = function(x,n){sample(x,n,replace = TRUE)}
n=5
m=500

y5= matrix(muestra(poblacion,n*m), ncol=n)


p5=apply(y5,1,sum)/n

par(mfrow = c(1, 2))
hist(p5,main = "Proporción, n=5",ylab = "Frecuencia", xlab = "p", col="#9fc3e7")
qqnorm(p5,main = " ", ylab = "Percenties Muestra", xlab = "Percentiles Normal", col="#0051a1")

S5=shapiro.test(p5)
S5


    Shapiro-Wilk normality test

data:  p5
W = 0.71826, p-value < 2.2e-16

#n=10

poblacion = c(rep(1,900),rep(0,100))
muestra = function(x,n){sample(x,n,replace = TRUE)}
n=10
m=500

y5= matrix(muestra(poblacion,n*m), ncol=n)


p6=apply(y5,1,sum)/n

par(mfrow = c(1, 2))
hist(p6,main = "Proporción, n=10",ylab = "Frecuencia", xlab = "p", col="#9fc3e7")
qqnorm(p6,main = " ", ylab = "Percenties Muestra", xlab = "Percentiles Normal", col="#0051a1")

S6=shapiro.test(p6)
S6


    Shapiro-Wilk normality test

data:  p6
W = 0.83662, p-value < 2.2e-16

#n=15

poblacion = c(rep(1,900),rep(0,100))
muestra = function(x,n){sample(x,n,replace = TRUE)}
n=15
m=500

y5= matrix(muestra(poblacion,n*m), ncol=n)


p7=apply(y5,1,sum)/n

par(mfrow = c(1, 2))
hist(p7,main = "Proporción, n=15",ylab = "Frecuencia", xlab = "p", col="#9fc3e7")
qqnorm(p7,main = " ", ylab = "Percenties Muestra", xlab = "Percentiles Normal", col="#0051a1")

S7=shapiro.test(p7)
S7


    Shapiro-Wilk normality test

data:  p7
W = 0.89881, p-value < 2.2e-16

#n=20

poblacion = c(rep(1,900),rep(0,100))
muestra = function(x,n){sample(x,n,replace = TRUE)}
n=20
m=500

y5= matrix(muestra(poblacion,n*m), ncol=n)


p8=apply(y5,1,sum)/n

par(mfrow = c(1, 2))
hist(p8,main = "Proporción, n=20",ylab = "Frecuencia", xlab = "p", col="#9fc3e7")
qqnorm(p8,main = " ", ylab = "Percenties Muestra", xlab = "Percentiles Normal", col="#0051a1")

S8=shapiro.test(p8)
S8


    Shapiro-Wilk normality test

data:  p8
W = 0.92479, p-value = 4.04e-15

#n=30


poblacion = c(rep(1,900),rep(0,100))
muestra = function(x,n){sample(x,n,replace = TRUE)}
n=30
m=500

y5= matrix(muestra(poblacion,n*m), ncol=n)


p9=apply(y5,1,sum)/n

par(mfrow = c(1, 2))
hist(p9,main = "Proporción, n=30",ylab = "Frecuencia", xlab = "p", col="#9fc3e7")
qqnorm(p9,main = " ", ylab = "Percenties Muestra", xlab = "Percentiles Normal", col="#0051a1")

S9=shapiro.test(p9)
S9


    Shapiro-Wilk normality test

data:  p9
W = 0.95016, p-value = 6.12e-12

#n=50

poblacion = c(rep(1,900),rep(0,100))
muestra = function(x,n){sample(x,n,replace = TRUE)}
n=50
m=500

y5= matrix(muestra(poblacion,n*m), ncol=n)


p10=apply(y5,1,sum)/n

par(mfrow = c(1, 2))
hist(p10,main = "Proporción, n=50",ylab = "Frecuencia", xlab = "p", col="#9fc3e7")
qqnorm(p10,main = " ", ylab = "Percenties Muestra", xlab = "Percentiles Normal", col="#0051a1")

S10=shapiro.test(p10)
S10


    Shapiro-Wilk normality test

data:  p10
W = 0.97158, p-value = 2.865e-08

#n=60

poblacion = c(rep(1,900),rep(0,100))
muestra = function(x,n){sample(x,n,replace = TRUE)}
n=60
m=500

y5= matrix(muestra(poblacion,n*m), ncol=n)


p11=apply(y5,1,sum)/n

par(mfrow = c(1, 2))
hist(p11,main = "Proporción, n=60",ylab = "Frecuencia", xlab = "p", col="#9fc3e7")
qqnorm(p11,main = " ", ylab = "Percenties Muestra", xlab = "Percentiles Normal", col="#0051a1")

S11=shapiro.test(p11)
S11


    Shapiro-Wilk normality test

data:  p11
W = 0.97303, p-value = 5.744e-08

#n=100

poblacion = c(rep(1,900),rep(0,100))
muestra = function(x,n){sample(x,n,replace = TRUE)}
n=100
m=500

y5= matrix(muestra(poblacion,n*m), ncol=n)


p12=apply(y5,1,sum)/n

par(mfrow = c(1, 2))
hist(p12,main = "Proporción, n=100",ylab = "Frecuencia", xlab = "p", col="#9fc3e7")
qqnorm(p12,main = " ", ylab = "Percenties Muestra", xlab = "Percentiles Normal", col="#0051a1")

S12=shapiro.test(p12)
S12


    Shapiro-Wilk normality test

data:  p12
W = 0.98189, p-value = 6.999e-06

#n=200

poblacion = c(rep(1,900),rep(0,100))
muestra = function(x,n){sample(x,n,replace = TRUE)}
n=200
m=500

y5= matrix(muestra(poblacion,n*m), ncol=n)


p13=apply(y5,1,sum)/n

par(mfrow = c(1, 2))
hist(p13,main = "Proporción, n=200",ylab = "Frecuencia", xlab = "p", col="#9fc3e7")
qqnorm(p13,main = " ", ylab = "Percenties Muestra", xlab = "Percentiles Normal", col="#0051a1")

S13=shapiro.test(p13)
S13


    Shapiro-Wilk normality test

data:  p13
W = 0.98567, p-value = 7.823e-05

#n=500

poblacion = c(rep(1,900),rep(0,100))
muestra = function(x,n){sample(x,n,replace = TRUE)}
n=500
m=500

y5= matrix(muestra(poblacion,n*m), ncol=n)


p14=apply(y5,1,sum)/n

par(mfrow = c(1, 2))
hist(p14,main = "Proporción, n=500",ylab = "Frecuencia", xlab = "p", col="#9fc3e7")
qqnorm(p14,main = " ", ylab = "Percenties Muestra", xlab = "Percentiles Normal", col="#0051a1")

S14=shapiro.test(p14)
S14


    Shapiro-Wilk normality test

data:  p14
W = 0.99414, p-value = 0.0516

#n=1000

poblacion = c(rep(1,900),rep(0,100))
muestra = function(x,n){sample(x,n,replace = TRUE)}
n=1000
m=500

y5= matrix(muestra(poblacion,n*m), ncol=n)


p15=apply(y5,1,sum)/n

par(mfrow = c(1, 2))
hist(p15,main = "Proporción, n=1000",ylab = "Frecuencia", xlab = "p", col="#9fc3e7")
qqnorm(p15,main = " ", ylab = "Percenties Muestra", xlab = "Percentiles Normal", col="#0051a1")

S15=shapiro.test(p15)
S15


    Shapiro-Wilk normality test

data:  p15
W = 0.99491, p-value = 0.09843

library(knitr)
library(kableExtra)
library(tidyverse)
# Crear la tabla
d <- data.frame(
  n = c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500, 1000),
  valor.p = c(S5$p.value, S6$p.value, S7$p.value, S8$p.value, S9$p.value, S10$p.value, S11$p.value, S12$p.value, S13$p.value, S14$p.value, S15$p.value),
  normalidad = c("no", "no", "no", "no", "no", "no", "no", "no", "no", "si","si")
)

# Crear la tabla con kable y kableExtra
styled_table <- d %>%
  mutate(
    normalidad = ifelse(normalidad == "si", "Sí", "No")
  ) %>%
  kable("html") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE) %>%
  column_spec(1, width = "70px") %>%
  column_spec(2, width = "100px") %>%
  column_spec(3, width = "100px")

styled_table

n	valor.p	normalidad
5	0.0000000	No
10	0.0000000	No
15	0.0000000	No
20	0.0000000	No
30	0.0000000	No
50	0.0000000	No
60	0.0000001	No
100	0.0000070	No
200	0.0000782	No
500	0.0515968	Sí
1000	0.0984271	Sí

En las graficas como en la tabla anterior se puede concluir que en la prueba de shapiro cuando el valor de p es mayor a 0.05 no se rechaza la hipotesis de normalidad, tambien en la tabla muestra que las proporciones para todos los n ninguno es normal menos el n=500 que si esta dentro de la normalidad, pero decidi agregarun n=1000 que si es normal, para esto se necesitan un mayor tamaño de muestra para encontrar la normalidad.

Problema 4

ESTIMACIÓN INTERVALOS DE CONFIANZA BOOTSTRAP

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∗, 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)

Construya el intervalo de confianza por los dos métodos y compare los resultados obtenidos. Comente los resultados. Confiaría en estas estimaciones?

7. Desarrollo

# Se desea estimar intervalos de confianza para la media de una muestra de datos.

# Datos de la muestra

x = c(7.69, 4.97, 4.56, 6.49, 4.34, 6.24, 4.45) # vector con los datos de la muestra.

# Número de muestras BOOTSTRAP

n_muestras = 1000

# Inicializar una matriz para almacenar las medias de BOOSTRAP.

medias_bootstrap = matrix(NA, nrow = n_muestras, ncol = length(x)) # x: valor de los datos de la muestra.

# Realizar el métosdo BOOTSTRAP

set.seed(418) # semilla para reproducibilidad.
for (i in 1:n_muestras) { # el "FOR" se utiliza para generar muestras, a partir de los datos de la muestra original que lleva el nombre de "X".
  muestra_bootstrap = sample(x, replace = TRUE) 
  medias_bootstrap[i,] = muestra_bootstrap # Es la matriz en donde se guarda la media de la muestra.
  
  
}

# Calcular las medias por filas, correspondiente a cada muestra Bootstrap generado.

medias_por_muestras <- apply(medias_bootstrap, 1, mean)

# Calcular intervalos de confianza para ambos métodos.

metodo_1 <- quantile(medias_por_muestras, probs = c(0.025, 0.975)) # el intervalo de confianza es calculado utilizando la función "quantile", con el cual se obtienen los cuartiles 2.5% y 97.5% de cada media por muestra.

metodo_1

    2.5%    97.5% 
4.737143 6.490143

metodo_2 <- c(2 * mean(medias_por_muestras) - metodo_1 [2], 2 * mean(medias_por_muestras) - metodo_1 [1]) # Se calculan los límites del intervalo. Se obtiene la media de las medias.

metodo_2

   97.5%     2.5% 
4.598109 6.351109

# Calcular intervalo de confianza para la media.

media_metodos = c(mean(medias_por_muestras) - sd(medias_por_muestras), mean(medias_por_muestras) + sd(medias_por_muestras))

# Graficar el histograma con líneas para intervalos de confianza

hist(medias_por_muestras, las = 1, main = "Histograma con intervalos de confianza y medida", xlab = "Media Bootstrap", ylab = "Frecuencia", col = "grey")

abline(v = metodo_1, col = "red", lwd = 2, lty = 2, label = "Método 1")
abline(v = metodo_2, col = "blue", lwd = 2, lty = 2, label = "Método 2")
abline(v = mean(medias_por_muestras), col = "green", lwd = 2, lty = 2, label = "Media")

legend("topright", legend = c("Método 1", "Método 2", "Media"), col = c("red", "blue", "green"), lty = 2, lwd = 2)

8. Nota

las muestras bootstrap se pueden obtener a partir de muestreo aleatorio con repetición (o tambien llamado con sustitución).
Entregable : enlace en RPubs con informe 4.
funciones recomendadas : sample() , apply(), quantile().
Problema tomado de Navidi(2006).

9. Análisis

Con el histograma anterior se puede ver la distribucion de las medias obtenidas basadas en la muestras bootstrap. La línea verde la cual se encuentra cerca al 5.5, representa la media de las medias obtenidas a través de las muestras bootstrap; por otro lado, las líneas roja y azul muestran los límites de los intervalos de confianza que fueron obtenidos a través de los dos métodos.

Teniendo en cuenta la ubicación de las líneas roja y azul, correspondiente al método 1 y al método 2 respectivamente, muestran que los intervalos de confianza para ambos métodos resultan ser similares y se encuentran ubicados al lado y lado de la media de las medias obtenidas a través de bootstrap.

Una vez comparados los dos métodos, ambos dan resultados positivos debido a que sus intervalos son razonables en cuanto a la estimación de la media poblacional.

Con una confianza del 95% se espera que la media de la eficiencia, se encuentre para el primer metodo entre 4,74 y 6,49. Para el segundo metodo estaria entre 4,60 y 6,35.

Actividad 2

Juan Camilo Gonzalez Victoria y Gabriel Orozco Castaño

2024-03-05

Problema 1

ESTIMACIÓN DE \(\pi\)

Conclusión:

Problema 2

PROPIEDADES DE LOS ESTIMADORES

4. Desarrollo

4.1. Muestra n = 20

4.2. Muestra n = 50

4.3. Muestra n = 100

4.4. Muestra n = 1000

4. Consistencia

4. Sesgo

4. Eficiencia

4.5. Conclusión - Problema 2

Problema 3

TEOREMA DEL LÍMITE CENTRAL

Problema 4

ESTIMACIÓN INTERVALOS DE CONFIANZA BOOTSTRAP

7. Desarrollo

8. Nota

9. Análisis