Problema 1

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 π.¿Cuántos de los puntos están dentro del círculo? ¿Cuál es su estimación de π?

iterations <- c(1000,10000,100000)
for (n in iterations) {
  x <- runif(n,min=0, max=1)
  y <- runif(n,min=0, max=1)
  inside_circle <- sum((x-0.5)^2+(y-0.5)^2 <0.25)
  PI_estimado =4*(inside_circle)/n
  
  cat("Numero de iteraciones N=", n, "\n")
  cat("Puntos dentro del circulo", inside_circle,"\n")
  cat("PI estimado:", PI_estimado, "\n")
  
}
## Numero de iteraciones N= 1000 
## Puntos dentro del circulo 794 
## PI estimado: 3.176 
## Numero de iteraciones N= 10000 
## Puntos dentro del circulo 7827 
## PI estimado: 3.1308 
## Numero de iteraciones N= 1e+05 
## Puntos dentro del circulo 78437 
## PI estimado: 3.13748

R// Como se observa en los resltados, conforme se aumenta el numero de muestras, nos acercamos al valor de π.

Problema 2

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.En cada caso evalue las propiedades de insesgadez, eficiencia y consistencia.

iterations <- c(20,50,100,1000)

theta=10

for(n in iterations){
    Theta1<-c()
    Theta2<-c()
    Theta3<-c()
    Theta4<-c()
  for (i in 1:n) {
    z <- rexp(4,rate=1/theta)
    theta1=((z[1]+z[2])/6) + ((z[3]+z[4])/3)
    theta2=((z[1]+2*z[2]+ 3*z[3]+4*z[4])/5)
    theta3=((z[1]+z[2]+z[3]+z[4])/4)
    theta4=(min(z)+max(z))/2
    Theta1 <-append(Theta1,theta1)
    Theta2 <-append(Theta2,theta2)
    Theta3 <-append(Theta3,theta3)
    Theta4 <-append(Theta4,theta4)
  }
  if  (n==20) {
    T20=data.frame(Theta1,Theta2,Theta3,Theta4)
  } else if (n==50) {
    T50=data.frame(Theta1,Theta2,Theta3,Theta4)
  } else if (n==100) {
    T100=data.frame(Theta1,Theta2,Theta3,Theta4)
  } else {
    T1000=data.frame(Theta1,Theta2,Theta3,Theta4)
  }
}

Boxplot 20

boxplot(T20, las=1, main="Comparacion estimadores con n=20")
abline(h=10, col="red")

apply(T20,2,mean) # Valores de las medias
##    Theta1    Theta2    Theta3    Theta4 
##  9.943013 20.010771 10.397246 12.423361
apply(T20,2,var) # Valores de la varianza
##   Theta1   Theta2   Theta3   Theta4 
## 19.70470 78.06729 20.73090 26.03873

R// En el caso de 20 iteraciones los estimadores con mejores resultados son θ3 y θ1 bastante cerca del parametro, se consideran insesgados por acercarse al parametro y presentar la menor varianza. θ2 es el parametro en este caso con mayor sesgo e ineficiencia.

Boxplot 50

boxplot(T50, las=1, main="Comparacion estimadores con n=50")
abline(h=10, col="red")

apply(T50,2,mean) # Valores de las medias
##    Theta1    Theta2    Theta3    Theta4 
##  9.111775 17.800551  9.377365 11.265525
apply(T50,2,var) # Valores de la varianza
##   Theta1   Theta2   Theta3   Theta4 
## 19.17926 76.30731 17.73417 36.28483

R// En el caso de 50 iteraciones el estimador con mejor resultado es θ3 al acercarse al parametro, se considera insesgado por ser el vañpr más cercanos al parametro y presentar la menor varianza. θ2 es el parametro en este caso con mayor sesgo e ineficiencia.

Boxplot 100

boxplot(T100, las=1, main="Comparacion estimadores con n=100")
abline(h=10, col="red")

apply(T100,2,mean) # Valores de las medias
##    Theta1    Theta2    Theta3    Theta4 
##  9.514471 19.081219  9.534484 11.173534
apply(T100,2,var) # Valores de la varianza
##   Theta1   Theta2   Theta3   Theta4 
## 23.37013 94.80073 21.19422 31.66902

R// En el caso de 100 iteraciones el estimador con mejor resultado son θ3 y θ1 bastante cerca del parametro y menor desviacion estandar, ambos se se consideran insesgados por ser los valores más cercanos al parametro y presentar la menor desviacion estandar (por ende varianza). θ2 es el parametro en este caso con mayor sesgo e ineficiencia.

Boxplot 1000

boxplot(T1000, las=1, main="Comparacion estimadores con n=1000")
abline(h=10, col="red")

apply(T1000,2,mean) # Valores de las medias
##    Theta1    Theta2    Theta3    Theta4 
##  9.949745 19.906519  9.925911 11.680674
apply(T1000,2,var) # Valores de la varianza
##    Theta1    Theta2    Theta3    Theta4 
##  27.85814 120.95821  24.68479  41.04717

R// En el caso de 1000 iteraciones el estimador con mejor resultado son θ3 y θ1 bastante cerca del parametro, ambos se consideran insesgados por ser los valores más cercanos al parametro y presentar la menor desviacion estandar (por ende varianza). θ2 es el parametro en este caso con mayor sesgo e ineficiencia.

R// Los parametros θ3 y θ1 parecen ser consistentes ya que parecen acercarse al valor real del parametro, 10. Con respecto θ4 se deben evaluar más casos con el fin de observar si presenta una aparente consistencia.y para θ2 no se observo consistencia ya que no se aproximo al valor del parametro con el aumento de muestras.

Problema 3

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:

Numerales a, b y c.

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

b. 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 500500 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

N <- rep(c(0,1),each=500) #poblacion inicial 50% enfermos

p_estimador <- function(size_sample){
  sum(sample(N,size =size_sample))/size_sample
}

caso70 <- numeric(500)

n=500

for (i in 1:n) {
  caso70[i] = p_estimador(70)
}

hist(caso70, main="Histograma escenario 70 veces", breaks=20, col="lightblue")

media <- mean(caso70)
varianza <- var(caso70)
valor_p <- shapiro.test(caso70)$p.value
cat("media ", media, "varianza ", varianza, "valor_p ", valor_p)
## media  0.4956571 varianza  0.003562141 valor_p  0.01264958

R// Como se observa parece buscar cierta simetria en el histograma, donde el estimador parece tender a estar en el rango de 0.45 y 0.55, con una media de 0.5 y un P-value mayor a 0.05, teniendo entonces una distribucion de tendencia normal.

Numeral d

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

lista_muestreos <- c(5,10,15,20,30,50,60,100,200,500)
 
for(i in lista_muestreos) {
  caso <- replicate(500, p_estimador(i))
  media <- mean(caso)
  varianza <- var(caso)
  valor_p <- shapiro.test(caso)$p.value
  cat("Para el caso de ", i," muestras tenemos una media de ", media," una varianza de ",varianza," un valor_p de ",valor_p,"\n")
  qqnorm(caso)
  qqline(caso)
  hist(caso, breaks=15)
}
## Para el caso de  5  muestras tenemos una media de  0.4932  una varianza de  0.052619  un valor_p de  1.649693e-14

## Para el caso de  10  muestras tenemos una media de  0.5074  una varianza de  0.02517559  un valor_p de  7.62894e-10

## Para el caso de  15  muestras tenemos una media de  0.4998667  una varianza de  0.01785792  un valor_p de  1.214309e-07

## Para el caso de  20  muestras tenemos una media de  0.4929  una varianza de  0.01188837  un valor_p de  4.280783e-06

## Para el caso de  30  muestras tenemos una media de  0.5008667  una varianza de  0.008364868  un valor_p de  0.0001380719

## Para el caso de  50  muestras tenemos una media de  0.49744  una varianza de  0.005340127  un valor_p de  0.002088048

## Para el caso de  60  muestras tenemos una media de  0.5012  una varianza de  0.003898579  un valor_p de  0.008482158

## Para el caso de  100  muestras tenemos una media de  0.49812  una varianza de  0.001990046  un valor_p de  0.05278609

## Para el caso de  200  muestras tenemos una media de  0.49942  una varianza de  0.001117699  un valor_p de  0.04974175

## Para el caso de  500  muestras tenemos una media de  0.498452  una varianza de  0.0002337913  un valor_p de  0.5959283

R// Se observa de manera grafica en los histogramas que conforme aumenta el valor de las muestras parece acercarse a un distribución normal, especialmente desde n=50. Este comportamiento tambien se puede observar en el grafico Q-Q normal, sin embargo visualmente este se acerca más al comportamiento normal cuando las muestras son iguales o mayores a 100. Mientras que observando el P_value se observa que conforme aumenta este numero de muestras este valor crece sin embargo para este ejercicio ninguno de ellos es mayor a 0.05

Numeral e

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.

Lote 10% plantas enfermas

N <- rep(c(0,1),c(900,100)) #poblacion inicial 10% enfermas

for(i in lista_muestreos) {
  caso <- replicate(500, p_estimador(i))
  media <- mean(caso)
  varianza <- var(caso)
  valor_p <- shapiro.test(caso)$p.value
  if (round(valor_p, digits = 2) > 0.05) {
    normalidad = "normal"
  } else {
    normalidad = "no normal"
  }
  cat("Para el caso de ", i," muestras tenemos una media de ", media," una varianza de ",varianza," un valor_p de ",valor_p,"indicando una distribucion",normalidad,"\n")
  qqnorm(caso)
  qqline(caso)
  hist(caso, breaks=15)
}
## Para el caso de  5  muestras tenemos una media de  0.09  una varianza de  0.01649299  un valor_p de  1.592752e-29 indicando una distribucion no normal

## Para el caso de  10  muestras tenemos una media de  0.0984  una varianza de  0.008213868  un valor_p de  6.992389e-22 indicando una distribucion no normal

## Para el caso de  15  muestras tenemos una media de  0.09946667  una varianza de  0.006154246  un valor_p de  3.737883e-18 indicando una distribucion no normal

## Para el caso de  20  muestras tenemos una media de  0.1025  una varianza de  0.004307365  un valor_p de  6.044643e-15 indicando una distribucion no normal

## Para el caso de  30  muestras tenemos una media de  0.0962  una varianza de  0.002450461  un valor_p de  4.384994e-12 indicando una distribucion no normal

## Para el caso de  50  muestras tenemos una media de  0.10368  una varianza de  0.001757974  un valor_p de  7.077848e-09 indicando una distribucion no normal

## Para el caso de  60  muestras tenemos una media de  0.0994  una varianza de  0.001497079  un valor_p de  4.329894e-07 indicando una distribucion no normal

## Para el caso de  100  muestras tenemos una media de  0.1  una varianza de  0.000907014  un valor_p de  8.242194e-05 indicando una distribucion no normal

## Para el caso de  200  muestras tenemos una media de  0.10011  una varianza de  0.0003552484  un valor_p de  0.0008013548 indicando una distribucion no normal

## Para el caso de  500  muestras tenemos una media de  0.09982  una varianza de  8.863287e-05  un valor_p de  0.01737234 indicando una distribucion no normal

Lote 90% plantas enfermas

N <- rep(c(0,1),c(100,900)) #poblacion inicial 90% enfermas

for(i in lista_muestreos) {
  caso <- replicate(500, p_estimador(i))
  media <- mean(caso)
  varianza <- var(caso)
  valor_p <- shapiro.test(caso)$p.value
  if (round(valor_p, digits = 2) > 0.05) {
    normalidad <- "normal"
  } else {
    normalidad <- "no normal"
  }
  cat("Para el caso de ", i," muestras tenemos una media de ", media," una varianza de ",varianza," un valor_p de ",valor_p,"indicando una distribucion",normalidad,"\n")
  qqnorm(caso)
  qqline(caso)
  hist(caso, breaks=15)
}
## Para el caso de  5  muestras tenemos una media de  0.902  una varianza de  0.0175511  un valor_p de  9.619702e-29 indicando una distribucion no normal

## Para el caso de  10  muestras tenemos una media de  0.8956  una varianza de  0.007876393  un valor_p de  4.50847e-21 indicando una distribucion no normal

## Para el caso de  15  muestras tenemos una media de  0.8941333  una varianza de  0.006120045  un valor_p de  4.733804e-17 indicando una distribucion no normal

## Para el caso de  20  muestras tenemos una media de  0.9015  una varianza de  0.004812375  un valor_p de  1.653628e-15 indicando una distribucion no normal

## Para el caso de  30  muestras tenemos una media de  0.8994  una varianza de  0.002905451  un valor_p de  1.059869e-11 indicando una distribucion no normal

## Para el caso de  50  muestras tenemos una media de  0.90108  una varianza de  0.001623681  un valor_p de  4.701037e-08 indicando una distribucion no normal

## Para el caso de  60  muestras tenemos una media de  0.8978  una varianza de  0.001139662  un valor_p de  1.239823e-08 indicando una distribucion no normal

## Para el caso de  100  muestras tenemos una media de  0.9006  una varianza de  0.0007010421  un valor_p de  3.244818e-06 indicando una distribucion no normal

## Para el caso de  200  muestras tenemos una media de  0.90018  una varianza de  0.0003875427  un valor_p de  0.001709697 indicando una distribucion no normal

## Para el caso de  500  muestras tenemos una media de  0.900612  una varianza de  9.965677e-05  un valor_p de  0.04164064 indicando una distribucion no normal

R// Para ambos casos 10% y 90% de las plantas están enfermas, de manera visual los histogramas parecen tener una tendencia simétrica solo en los casos por encima de 100 muestras (100, 200, y 500), lo cual coincide con lo visualizado en el gráfico Q-Q normal específicamente para el caso de 500. Sin embargo para el P-value no se observa en ninguno de los casos valores mayores 0.05 por lo cual no siguen tendencias normales, aunque si se observa que conforme aumenta el numero de muestras el p-value incrementa conforme lo hace el número de muestreos.

Problema 4

Es necesario extraer un gran número de muestras (suponga k = 1000). Para cada una de las muestra bootstrap obtenidas se calcula la media, obteniéndose un valor para cada muestra. El intervalo de confianza queda conformado por los percentiles P2.5 y P97.5.

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

muestra =c(7.69, 4.97, 4.56, 6.49, 4.34,6.24,4.25)

muestra_b=sample(muestra, 7000, replace=TRUE)

matrix_boot= matrix(muestra_b, nrow=1000, ncol=7)
medias=rowMeans(matrix_boot)

cat("Media ", mean(medias), " varianza", var(medias), " p_value ",shapiro.test(medias)$p.value)
## Media  5.528817  varianza 0.2242541  p_value  0.005931473
intervalo1 = quantile(medias, probs = c(0.025,0.975))
intervalo2=c(2*mean(medias)-quantile(medias,0.975),2*mean(medias)-quantile(medias,0.025))
print(intervalo1) 
##     2.5%    97.5% 
## 4.679536 6.498821
print(intervalo2)
##    97.5%     2.5% 
## 4.558813 6.378099
hist(medias, las=1, main="Histograma y estimado intervalos confianza", col="green", breaks=20)
abline(v=intervalo1, col="red",lwd=2)
abline(v=intervalo2, col="blue",lwd=2)

qqnorm(medias)
qqline(medias)

R// Ambos resultados son bastante cercanos en el rango entre los percentiles 2.5% y 97.5%. De igual forma ambos percentiles se encuentran entre el maximo (7.69) y minimo del vector inicial (4.25), la forma de su histograma parece seguir una tendencia normal al igual que el grafico Q-Q normal. El unico parametro que invalida la hipotesis es el P-value inferior a 0.05. De manera preliminar parecen un buen estimado de este intervalo de confianza, sin embargo se recomendaría un mayor numero de muestreos con el objetivo de buscar un P-value cercano o mayor a 0.05.