Estimación del valor de π.
La siguiente figura sugiere como estimar el valor de con una simulación. En la figura, un círculo 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 π.
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).
Genere 1000 coordenadas y: Y1,…, Yn, utilizando nuevamente la distribución uniforme con valor mínimo de 0 y valor máximo de 1.
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.
¿Cuántos de los puntos están dentro del círculo? ¿Cuál es su estimación de π?
## [1] 0.2875775 0.7883051 0.4089769 0.8830174 0.9404673 0.0455565 0.5281055
## [8] 0.8924190 0.5514350 0.4566147
## [1] 0.27362273 0.59386693 0.16018481 0.85343024 0.84773916 0.47788681
## [7] 0.77369212 0.29540008 0.06562811 0.44053312
## [1] 0.09636998 0.09193085 0.12375956 0.27161527 0.31493395 0.20700789
## [7] 0.07569730 0.19585383 0.19132450 0.00541859
## [1] 800
## [1] 3.2
De acuerdo con lo planteado en los pasos del ejercicio, se generaron 1000 coordenadas utilizando un valor mínimo de 0 y un máximo de 1. Se pueden observar las 10 primeras coordenadas del eje x relacionadas con las 10 primeras coordenadas del eje y. Una vez establecido esto, se procedió a calcular las distancias de las coordenadas a través de (Xi−0.5)2 + (Yi−0.5)2, y se muestran los 10 primeros resultados redondeados a 10 cífras decimales. Realizado esto, se establece la condición numérica y se calculan los valores inferiores a 0.25, estos son los valores que se encuentran dentro del área del círculo. Finalmente, se procede a la estimación del valor de π de acuerdo con la muestra de 1000 puntos.
A partir de la simulación efectuada por medio de 1000 coordenadas, se observa lo siguiente: el número de puntos dentro del círculo es equivalente a 800, mientras el número de puntos fuera del círculo es 200. Para 1000 coordenadas el valor estimado de π es 3,2. Para una muestra de 10000 coordenadas se obtienen 7894 puntos dentro del círculo, y 2106 puntos fuera del círculo, además de una estimación del valor de π equivalente a 3.1576, mucho más cercana al valor irracional del número. Para 1000000 de coordenadas la estimación de π es 3.141408. En términos generales, entre más coordenadas se tomen mayor será la aproximación al valor exacto de π, o entre mayor sea la muestra más cerca se estará del valor.
n = 1000; 10000; 1000000
x = runif(n)
y = runif(n)
distancia = (x - 0.5)^2 + (y - 0.5)^2
t = as.numeric(distancia < 0.25)
sum(t)
pi = sum((t)/n)*4
fx = function(x,y){sqrt((x-0.5)^2 + (y-0.5)^2)}
fx = function(row) {distancia = sqrt(row[‘x’]^2 + row[‘y’]^2) return(distancia) }
Eficiencia, consistencia e insesgadez de estimadores.
Sean X1, X2, X3 y X4, una muestra aleatoria de tamaño n = 4 cuya población la conforma una distribución exponencial con parámetro tetha desconocido. Determine las características de cada uno de los siguientes estimadores propuestos.
Genere una muestras de n = 20, 50, 100 y 1000 para cada uno de los estimadores planteados.
En cada caso evalue las propiedades de insesgadez, eficiencia y consistencia
Suponga un valor para el parámetro θ
θ1ˆ=X1+X2/6 + X3+X4/3
θ2ˆ=(X1+2X2+3X3+4X4)5
θ3ˆ=X1+X2+X3+X4/4$
θ4ˆ=min{X1,X2,X3,X4}+max{X1,X2,X3,X4}2
set.seed(123)
n = 100
x1 = rexp(n, 0.5)
x2 = rexp(n, 0.5)
x3 = rexp(n, 0.5)
x4 = rexp(n, 0.5)
base = data.frame(x1, x2, x3, x4)
fxa = function(x) {(x[1]+x[2])/6+(x[3]+x[4])/3}
fxb = function(x) {(x[1]+2*x[2]+3*x[3]+4*x[4])/5}
fxc = function(x) (x[1]+x[2]+x[3]+x[4])/4
fxd = function(x) {return(min(x)+max(x)/2)}
T1 = apply(base,1,fxa)
T2 = apply(base,1,fxb)
T3 = apply(base,1,fxc)
T4 = apply(base,1,fxd)
A = data.frame(T1,T2,T3,T4)
boxplot(A, main = "Comparación de Estimadores", ylab = "Valores", col = "lightblue")
abline(h = 2, col = 'red')## T1 T2 T3 T4
## Min. :0.1481 Min. :0.2794 Min. :0.1339 Min. :0.1943
## 1st Qu.:1.3125 1st Qu.:2.5962 1st Qu.:1.3827 1st Qu.:1.6547
## Median :1.8791 Median :3.6031 Median :1.9298 Median :2.3410
## Mean :1.9767 Mean :3.9133 Mean :1.9861 Mean :2.5437
## 3rd Qu.:2.6029 3rd Qu.:4.8720 3rd Qu.:2.4674 3rd Qu.:3.1224
## Max. :4.6288 Max. :9.9914 Max. :4.9674 Max. :7.4779
## T1 T2 T3 T4
## T1 0.8288984 1.717970 0.7517189 0.9863185
## T2 1.7179699 3.756859 1.4921595 1.9490801
## T3 0.7517189 1.492159 0.7913008 1.0380056
## T4 0.9863185 1.949080 1.0380056 1.5916324
Para este ejercicio se asignó el valor de 2 al parámetro θ y 0.5 al valor de λ. Se asociaron 20, 50, 100 y 1000 números expoenciales a cada selección aleatoria, respectivamente. Una vez hecho esto, se crearon cuatro funciones asociadas con cuatro vectores para cada estimador proupuesto. Hecho esto, se inició el calculo de los estimadores a través de la función apply. Ya calculados los estimadores se comenzó la graficación de la información por medio de un diagrama de Cajas y Bigotes. Para la comprensión de los datos se corrieron las funciones summary y var. Este procedimiento se repitió cuatros veces con el fin de evaluar los resulados de cada número exponencial.
De acuerdo con el objetivo del ejercicio, se identificaron las propiedades de insesgadez, eficiencia y consistencia.Para una muestra de 100 números explenciales se generaron 4 estimadores: T1, T2, T3 y T4 con un parámetro θ = 2. El estimador que presentó mejor rendimiento fue T3, seguido de T1, T4 y T2. El estimador T3 presentó un grado de insesgadez muy alto ya que obtuvo una media de 1.9861. En cuanto a la eficiencia, T3 presentó una varianza de 0.7913 lo que establece un buen rendimiento. Finalmente, el estimador T3 se asoció con el valor de consistencia debido a que presentó mejor rendimiento al momento de incrementar el valor exponencial de la muestra. En términos generales, los estimadores presentaron mejores rendimientos una vez que se incrementaba el valor de la muestra.
Al comparar los dos gráficos se observa como incrementa el rendimiento de los estimadores una vez la muestra aumenta.
library(e1071)
library(nortsTest)
Teorema del Límite Central.
Esta técnica es una 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.
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. 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.
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)
muestra_proporcion = function(n) {
muestra = sample(c(0, 1), size = n, replace = TRUE, prob = c(0.5, 0.5))
prop_muestra = mean(muestra)
return(prop_muestra)
}
n <- c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)
estimadores = list()
pvalores = c()
for (i in n) {
estimadores[[as.character(i)]] <- replicate(500, muestra_proporcion(i))
prueba <- shapiro.test(estimadores[[as.character(i)]])
pvalores <- c(pvalores, prueba$p.value)
}
boxplot(estimadores, main = "Estimadores por tamaño de muestra",
xlab = "Tamaño de muestra", ylab = "Estimador", varwidth = TRUE, col = "lightblue")plot(n, pvalores, main = "P-valores por tamaño de muestra",
xlab = "Tamaño de muestra", ylab = "P-valor", pch = 19, col = "red")par(mfrow = c(2, 5))
for (i in n) {
qqnorm(estimadores[[as.character(i)]], main = paste("Gráfico QQ normal para n =", i))
}Para la realización de una simulación que genere una población n = 1000 donde el porcentaje de individuos enfermos corresponda al 50% se establece la siguiente función:
muestra_proporcion = function(n) {muestra = sample (c (0,1), size = n, replace = TRUE, prob = c (0.5, 0.5)) prop_muestra = mean(muestra) return(prop_muestra}
Una vez definido esto, se crea una lista vacía para almacenar los estimadores, además de crear un vector vacío para los p-valores originados por la prueba de normalidad. Por otra parte, se crea un vector n donde se registran los valores pedidos en el ejercicio: 5, 10, 15, 20, 30, 50, 60, 100, 200 y 500. Finalmente, se obtienen los estimadores y se calcula el p-valor. Por último, se visualizó la información por medio de gráficos de Cajas y Bigotes, Dispersión y QQplot.
El anterior procedimiento se repitió al momento de ejecutar el punto e. del ejercicio, esto es, la simulación de lotes con el 10% de plantas enfermas, y 90% de plantas enfermas. En este caso se uso la función base:
muestra_proporcion = function(n) {muestra = sample (c (0,1), size = n, replace = TRUE, prob = c (0.1, 0.9)) prop_muestra = mean(muestra) return(prop_muestra}
muestra_proporcion = function(n) {muestra = sample (c (0,1), size = n, replace = TRUE, prob = c (0.9, 0.1)) prop_muestra = mean(muestra) return(prop_muestra}
Se puede ver como el comportamiento de los promedios se aproxima mejor a una distribución normal cuando incrementa el tamaño de la muestra, sin importar la forma de la distribución original. De acuerdo con lo solicitado, la distribución original es de carácter binomial por lo que pudo presentar un comportamiento simétrico o asimétrico el cual dependió de la probabilidad establecida de 0.5, 0.1 y 0.9. Al momento de evaluar la información registrada en los diagramas, se observa que el boxplot indica que la variabilidad de los estimadores disminuye cuando incremental la muestra. Al revisar los diagramas de dispersión, se ve como los p-valores incrementan en la medida que lo hace la muestra, lo que avala la hipótesis nula frente a la distribución normal de los estimadores en muestras altas, 500, por ejemplo, además de ser simétricos y eficientes (con menor varianza). Por último, el QQplot muestra como los puntos se alinean en la recta diagonal cuando el tamaño de la muestra es mayor.
El comportamiento de los estimadores puede verse claramente al comparar los tres casos probabilísticos planteados en el ejercicio 3.
Los siguientes histogramas muestran las distribuciones de los estimadores para replicate = 500, de acuerdo con los tres casos de probabilidad establecidos en el ejercicio.
## El porcentaje de plantas enfermas es 50.7 %
## [1] 0.507
## El porcentaje de plantas enfermas es 90.8 %
## [1] 0.908
## El porcentaje de plantas enfermas es 9.2 %
## [1] 0.092
library(magrittr)
library(ggpubr)
library(ggplot2)
library(stats)
library(moments)
library(e1071)
library(nortsTest)
Método Bootstrap.
A partir de los siguientes datos: 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 , regresándolo 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
Construya el intervalo de confianza por los dos métodos y compare los resultados obtenidos. Comente los resultados. ¿Confiaría en estas estimaciones?
Para la solución del ejercicio 4 use los dos métodos propuestos.
Primer método.
set.seed(123)
k = 1000
muestra = c(7.69, 4.97, 4.56, 6.49, 4.34, 6.24, 4.45)
medias_bootstrap = numeric(k)
for (i in 1:k) {
muestra_bootstrap = sample(muestra, replace = TRUE)
medias_bootstrap[i] = mean(muestra_bootstrap)
}
percentil_2.5 = quantile(medias_bootstrap, 0.025)
percentil_97.5 = quantile(medias_bootstrap, 0.975)
intervalo_confianza = c(percentil_2.5, percentil_97.5)
cat("Intervalo de confianza del 95%", intervalo_confianza[1], "-", intervalo_confianza[2], "\n")## Intervalo de confianza del 95% 4.748393 - 6.508643
Segundo Método.
set.seed(123)
x = c(7.69, 4.97, 4.56, 6.49, 4.34, 6.24, 4.45)
boot = sample(x, 7000, replace = TRUE)
b = matrix(boot, nrow = 1000, ncol = 7)
mx = apply(b, 1, mean)
ic1 = quantile(mx, probs = c(0.025, 0.975))
print(ic1)## 2.5% 97.5%
## 4.721429 6.418643
## 97.5% 2.5%
## 4.639526 6.336740
hist(mx, las = 1, main = '', ylab = '', col = 'lightblue' )
abline(v = ic1, col = 'red', lwd = 2)
abline(v = ic2, col = 'blue', lwd = 2)Para la solución del ejercicio 4 se creó un vector para almacenar las muestras boostrap con el fin de generar una muestra que permitiera calcular su valor medio. Una vez desarrallado esto, se procedió al cálculo del intervalo de confianza del 95%. Los resultados encontrados indican que existe un 95% de certeza de que los datos se encuentran dentro de un intervalo que inicia en 4.748393 y termina en 6.508643 (esto a través del primer método). Por otra parte, el segundo método indica que el 95% de certeza se encuentra entre 4.639526 y 6.336740. Estos resultados son muy similares lo que lleva a pensar que las estimación son correctas, además de apoyarse en la distribución normal de los estimadores (ver el histograma). En términos generales, la ubicación del parámetro real dentro del percentil 2.5% y 97.5% es muy similar en los 2 métodos.
Se realizó un análisis Bootstrap con el fin de establecer un intervalo de confianza del 95% capaz de estimar el grado de incertidumbre de la eficiencia de un combustible en una muestra aleatoria de camiones. Para esto, se generaron 1000 muestras a las que se calculó las medias para la construcción de un intervalo por medio dos métodos distintos. Se observa que la muestra elegida (y consignada en este informe) resultó ser representativa de la población ya que se ajustó de manera muy similar a los resultados alcanzados a partir de los dos métodos de cálculo. A partir del uso de técnicas de probabilidad, se ejecutaron diversas muestras con el fin de evaluar el comportamiento aleatorio del remuestreo con reemplazo. No obstante, los resultados fueron similares en todos los casos debido al tamaño de la muestra, 1000.
boot = sample(x, 7000, replace = TRUE)
b = matrix(boot, nrow = 1000, ncol = 7)