markdown
Este proyecto tiene como objetivo principal determinar diferentes intervalos de confianza que generar diferentes métodos para la estimación del parámetro \(\beta\) de la distribución Rayleigh. Podemos comparar estos métodos a través de simulaciones de Monte Carlo considerando tanto aproximaciones asintóticas como métodos basados en propiedades exactas de la distribución. Este tipo de ejercicios resulta fundamental en estadística matemática, ya que permite realizar un contraste entre teoría y práctica, al mismo tiempo que se determina la efectividad de distintos estimadores frente a variaciones en el tamaño y condiciones de la muestra.
La distribución Rayleigh\((\beta)\), tal que \(\beta > 0\), tiene aplicaciones en diversos contextos de modelamiento estadístico, por ejemplo en el análisis de señales, vibraciones y propagación de ondas. La función de densidad de la distribución Rayleigh(\(\beta\)) está dada por
\[ f_\beta(x) = \frac{x}{\beta^2}\, e^{-x^2 / 2\beta^2}, \qquad x > 0, \ \beta > 0. \].
En adición, la media y varianza para esta distribución son
\[ \mathbb{E}[X] = \beta \sqrt{\frac{\pi}{2}}, \qquad \operatorname{Var}(X) = \beta^2 \left( 2 - \frac{\pi}{2} \right) \],
respectivamente.
De manera más general, para \(j \geq 1\) y \(X \sim \text{Rayleigh}(\beta)\), tenemos que los momentos estan dados por
\[ \mathbb{E}[X^j] = \beta^j 2^{j/2}\,\Gamma\!\left(1 + \frac{j}{2}\right), \]
donde \(\Gamma(\cdot)\) representa la función Gamma. Más adelante, vamos a aprovechar los momentos de la distribución para la construccion de intervalos de confianza que nos permitirán hacer estimaciones del parámetro \(\beta\).
Una de las ventajas de la distribución Rayleigh es que puede simularse fácilmente a partir de variables aleatorias uniformes. De este modo, vamos a probar que si \(U \sim \text{Unif}(0,1)\), entonces
\[ X = \sqrt{-2\beta^2 \ln(1-U)} \sim \text{Rayleigh}(\beta). \]
Con base a esto, podemos generar una muestra de tamaño \(n\) de variables con la distribución de Rayleigh simulando las variables aleatorias \(U_1,\dots,U_n\) independientes e idénticamente distribuidas con la distribución uniforme en \((0,1)\) y aplicando la transformación anterior. Este proceso nos permitirá que como resultado obtengamos una muestra con la distribución deseada.
Para la verificación de la validez de la muestra generada con la distribución de Rayleigh se pueden emplear gráficos de cuantiles tambien conocidos como q-q plots, que permiten comparar los cuantiles de dos distribuciones diferentes. Para este caso, los q-q plots nos van a permitir comparar los cuantiles teóricos de la distribución Rayleigh con cuantiles obtenidos de una muestra generada. Para el análisis de estos gráficos, si los puntos se alinean cerca de la diagonal \(y = x\), podemos concluir que el algoritmo de generación produce datos consistentes con respecto a la distribución teórica. Esta comparación es importante ya que de obtener resultados positivos, podemos iniciar con los cálculos para construir los intervalos de confianza.
El interés principal del proyecto es estudiar la efectividad de diferentes procedimientos a la hora de estimar los intervalos de confianza para \(\beta\). Entre los métodos considerados se encuentran:
Intervalos basados en la media muestral y el Teorema del Límite Central.
Intervalos a partir de los cuadrados de las variables aleatorias, \(V_i = X_i^2\).
Intervalos basados en la mediana muestral, la distribución exacta de los estadísticos de orden y la distribución beta.
Procedimientos alternativos que pueden ser asintóticos o exactos.
Cada uno de estos métodos presenta ventajas y limitaciones, principalmente en la probabilidad estimada de contener el parámetro, el ancho de los intervalos y la incertidumbre en los estimadores producidos.
Para comparar los diferentes métodos vamos a recurrir a simulaciones de Monte Carlo con diferentes variaciones en las condiciones. Esta simulación consiste en repetir un experimento computacional una gran cantidad de veces, para este caso \(B=50{,}000\), \(n=50,100.200,500\) y \(\alpha=0.02\); con esto, la idea es calcular en cada repetición los intervalos de confianza. Con base en estas repeticiones se puede estimar:
La cobertura, es decir, la probabilidad estimada de que los intervalos contengan al verdadero valor de \(\beta\).
La longitud de los intervalos.
La incertidumbre asociada a estas estimaciones.
En resumen, con el proyecto buscamos justificar de manera teórica y computacional los métodos de generación de muestras para la distribución de Rayleigh. Además de determinar diversos métodos de construcción de intervalos de confianza para \(\beta\). Por otra parte, queremos implementar simulaciones Monte Carlo que permitan comparar estos métodos en términos de cobertura, longitud, incertidumbre y complejidad computacional. Finalmente, de estos resultados buscamos extraer conclusiones sobre la efectividad de cada método en función del tamaño muestral y los parámetros dados.
Primero, justifiquemos la validez del Algoritmo 1.
Calculemos \(F(x)\) para \(X \sim \text{Rayleigh}\).
Tomando \[ u = -\frac{t^{2}}{2\beta^{2}}, \quad \text{entonces } du = -\frac{t}{\beta^{2}}\, dt, \] luego
\[ F(x) = P(X \leq x) = \int_{0}^{x} \frac{t}{\beta^{2}} e^{-\tfrac{t^{2}}{2\beta^{2}}} \, dt. \]
Con el cambio de variable:
\[ F(x) = - \int_{0}^{-\tfrac{x^{2}}{2\beta^{2}}} e^{u} \, du = - e^{u} \Big|_{0}^{-\tfrac{x^{2}}{2\beta^{2}}}. \]
\[ F(x) = - \left[ e^{-\tfrac{x^{2}}{2\beta^{2}}} - 1 \right] = 1 - e^{-\tfrac{x^{2}}{2\beta^{2}}}. \]
Sean \(\; U_1, U_2, \ldots, U_n \;\; i.i.d. \;\; \text{Unif}(0,1)\). Entonces, tomemos
\[ X_i = F^{-1}(U_i) \]
De este modo,
\[ u_i = 1 - e^{-\tfrac{x_i^2}{2\beta^2}} \]
\[ e^{-\tfrac{x_i^2}{2\beta^2}} = 1 - u \]
\[ -\tfrac{x_i^2}{2\beta^2} = \ln(1-u) \]
\[ x_i^2 = -2\beta^2 \ln(1-u) \]
\[ x_i = \sqrt{-2\beta^2 \ln(1-u)}. \]
De esta manera,
\[ P(X_i \leq x) = P(F^{-1}(U_i) \leq x) = P(U_i \leq F(x)) = F(x), \]
ya que ; \(U_i \sim \text{Unif}(0,1)\).
De este modo, \(X_i \sim \text{Rayleigh}\) y \(X_1, \ldots, X_n\) es una muestra i.i.d. de la Rayleigh.
Para validar experimentalmente este algoritmo podemos construir un “q-q plot” que compara gráficamente los cuantiles de dos distribuciones. De este modo, si dos distribuciones son similares, la gráfica del “q-q plot” debe ser similar a la recta \(y=x\).
#Generamos los datos para U con distribución Unif(0,1)
set.seed(123)
u <- runif(2000)
Ahora, utilizando el algoritmo y tomando \(\beta = 2.5\), construimos una muestra de 2000 datos de la distribución Rayleigh(\(\beta = 2.5\)). En adición, para no imprimir todos los datos, podemos visualizar los primeros datos.
#Calculamos x con el algortimo.
x <- sqrt(-2 * (2.5^2) * log(1 - u))
head(x)
## [1] 2.0587744 4.4054077 2.5639330 5.1789600 5.9384649 0.7634366
Sabemos que el cuantil de \(p\) es \(\zeta_p \in Rango(x)\), tal que \(F(\zeta_p)=p\) .
Por lo tanto, \(F^{-1}(p) = \zeta_p\); es decir el cuantil teórico de \(p\) para la distribución Rayleigh con parámetro \(\beta = 2.5\) se define como
\[ \zeta_p = \sqrt{-2 \beta^2 \ln(1-p)}. \]
Con esto, calculando los cuantiles teóricos de probabilidad 0.01 (0.01) 0.99 de la Rayleigh(2.5), obtenemos
#Le asignamos a p los valores entre 0.01 y 0.99 con pasos de 0.01
p <- seq(0.01, 0.99, by = 0.01)
#Calculamos los cuantiles teóricos de p
q_teoricos <- sqrt(-2 *(2.5^2)* log(1 - p))
q_teoricos
## [1] 0.3544421 0.5025275 0.6170414 0.7143353 0.8007285 0.8794558 0.9524356
## [8] 1.0209163 1.0857640 1.1476090 1.2069270 1.2640875 1.3193846 1.3730572
## [15] 1.4253023 1.4762850 1.5261454 1.5750037 1.6229642 1.6701181 1.7165457
## [22] 1.7623186 1.8075009 1.8521503 1.8963190 1.9400551 1.9834022 2.0264010
## [29] 2.0690889 2.1115011 2.1536704 2.1956277 2.2374024 2.2790224 2.3205143
## [36] 2.3619036 2.4032152 2.4444728 2.4856999 2.5269191 2.5681529 2.6094232
## [43] 2.6507521 2.6921611 2.7336720 2.7753066 2.8170869 2.8590349 2.9011734
## [50] 2.9435251 2.9861135 3.0289626 3.0720974 3.1155433 3.1593269 3.2034758
## [57] 3.2480188 3.2929860 3.3384093 3.3843218 3.4307589 3.4777579 3.5253586
## [64] 3.5736033 3.6225373 3.6722092 3.7226715 3.7739805 3.8261975 3.8793891
## [71] 3.9336280 3.9889937 4.0455737 4.1034645 4.1627731 4.2236186 4.2861346
## [78] 4.3504709 4.4167971 4.4853064 4.5562199 4.6297927 4.7063213 4.7861538
## [85] 4.8697022 4.9574601 5.0500258 5.1481350 5.2527075 5.3649151 5.4862847
## [92] 5.6188618 5.7654792 5.9302305 6.1193671 6.3431812 6.6205720 6.9928741
## [99] 7.5871356
# Calculemos los cuantiles muestrales de la muestra x
q_muestral <- as.numeric(quantile(x, probs = seq(0.01,0.99,0.01)))
q_muestral
## [1] 0.3624833 0.5034645 0.6130286 0.7176474 0.7813318 0.8677542 0.9614868
## [8] 1.0046130 1.0719900 1.1434632 1.2056035 1.2574490 1.3191824 1.3642497
## [15] 1.4095320 1.4548969 1.5007856 1.5627778 1.6296334 1.6872768 1.7458009
## [22] 1.7789422 1.8277226 1.8679943 1.9121903 1.9612514 2.0002179 2.0367772
## [29] 2.0652507 2.1222085 2.1671050 2.2022726 2.2482956 2.3039340 2.3487761
## [36] 2.3955175 2.4307749 2.4703218 2.5022006 2.5453304 2.5808432 2.6061241
## [43] 2.6446752 2.6907335 2.7359755 2.7775615 2.8024640 2.8494769 2.8861122
## [50] 2.9220798 2.9677429 3.0227274 3.0634980 3.1099813 3.1492260 3.1853023
## [57] 3.2214730 3.2574490 3.3059246 3.3501356 3.3818135 3.4437398 3.4864606
## [64] 3.5245316 3.5837491 3.6235909 3.6690041 3.7226923 3.7880390 3.8420455
## [71] 3.9119919 3.9703862 4.0318531 4.0780043 4.1401344 4.1936708 4.2536101
## [78] 4.3304241 4.3887697 4.4714510 4.5272605 4.5741766 4.6334827 4.7205330
## [85] 4.8190363 4.8938393 4.9734095 5.0488152 5.1608983 5.2361698 5.3381034
## [92] 5.5469144 5.7821317 5.9132743 6.1232498 6.3098784 6.6157807 7.1556401
## [99] 7.6952652
#Vamos a generar el q-q plot.
q_teoricos <- as.numeric(q_teoricos); names(q_teoricos) <- NULL
q_muestral <- as.numeric(q_muestral); names(q_muestral) <- NULL
qqplot(q_teoricos, q_muestral, xlab="Cuantiles teóricos", ylab="Cuantiles muestrales")
abline(0,1,col="red")
De la gráfica, podemos observar que la recta que se forma con los puntos es muy similar a la recta \(y=x\). Esto quiere decir que los cuantiles dados por la muestra son casi iguales (o iguales) a los cuantiles teóricos calculados. Por lo tanto, la muestra tiene distribución muy cercana a la distribución \(Rayleigh(2.5)\)
Para n=500, tenemos que
#primero simulamos la muestra con la distribución uniforme.
set.seed(123)
u_500 <- runif(500)
#Calculamos x con el algortimo.
x_500 <- sqrt(-2 * (2.5^2) * log(1 - u_500))
#Luego, calculamos los cuantiles muestrales.
q_muestral_500 <- as.numeric(quantile(x_500, probs = seq(0.01,0.99,0.01)))
#Finalmente, generamos el q-q plot.
q_teoricos <- as.numeric(q_teoricos); names(q_teoricos) <- NULL
q_muestral <- as.numeric(q_muestral_500); names(q_muestral_500) <- NULL
qqplot(q_teoricos, q_muestral_500, xlab="Cuantiles teóricos", ylab="Cuantiles muestrales")
abline(0,1,col="red")
Para n=20000, tenemos que
#primero simulamos la muestra con la distribución uniforme.
set.seed(123)
u_20000 <- runif(20000)
#Calculamos x con el algortimo.
x_20000 <- sqrt(-2 * (2.5^2) * log(1 - u_20000))
#Luego, calculamos los cuantiles muestrales.
q_muestral_20000 <- as.numeric(quantile(x_20000, probs = seq(0.01,0.99,0.01)))
#Finalmente, generamos el q-q plot.
q_teoricos <- as.numeric(q_teoricos); names(q_teoricos) <- NULL
q_muestral <- as.numeric(q_muestral_20000); names(q_muestral_20000) <- NULL
qqplot(q_teoricos, q_muestral_20000, xlab="Cuantiles teóricos", ylab="Cuantiles muestrales")
abline(0,1,col="red")
A partir de estas gráficas, podemos evidenciar que para n=500, n=2000 y n=20000; las gráficas se comportan de manera muy similar. Sin embargo, a medida que aumenta la cantidad de datos, los puntos forman una recta más definida y similar a la recta \(y=x\). De este modo, cada gráfica muestra la cercanía de la distribución muestral con la distribución de \(Rayleigh(2.5)\), tal que se hace aún más similar al aumentar la cantidad de datos en la muestra.
Veamos que \[ \frac{\sqrt{n}}{S}\left(\overline{X} - \beta \sqrt{\frac{\pi}{2}}\right) \approx N(0,1). \]
Para \(X_i \sim \text{Rayleigh}(\beta)\), tenemos que \[ \mathbb{E}(X_i) = \beta \sqrt{\frac{\pi}{2}} = \mu. \]
Además, sea \[ S_n = X_1 + X_2 + \cdots + X_n. \]
Luego, \[ \frac{\sqrt{n}}{S}\left(\overline{X} - \beta \sqrt{\frac{\pi}{2}}\right) = \frac{\sqrt{n}}{S}\left(\frac{X_1 + X_2 + \cdots + X_n}{n} - \mu\right), \]
\[ = \frac{\sqrt{n}}{S}\left(\frac{X_1 + X_2 + \cdots + X_n - n\mu}{n}\right), \]
\[ = \frac{X_1 + X_2 + \cdots + X_n - n\mu}{\sqrt{n}S} = \frac{S_n - n\mu}{\sqrt{n}S}. \]
Utilizando el Teorema del Límite Central, sabemos que \[ \frac{S_n - n\mu}{\sqrt{n}S} \;\;\xrightarrow[n\to\infty]{}\;\; N(0,1). \]
De este modo, \[ \frac{S_n - n\mu}{\sqrt{n}S} = \frac{\sqrt{n}}{S}\left(\overline{X} - \beta \sqrt{\frac{\pi}{2}}\right) \approx N(0,1). \,\,\,\,\,\,\, (1) \]
Ahora, deduzcamos un intervalo de confianza para \(\beta\).
Sea \(Z \sim N(0,1)\), entonces, un intervalo del \(1-\alpha\) de confianza para \(\beta\) es
\[ P(-{\zeta}_{1-{\alpha}/2} \leq Z \leq {\zeta}_{1-{\alpha}/2}) = 1-{\alpha} \]
Teniendo en cuenta lo probado en (1), tenemos que
\[ P(-{\zeta}_{1-{\alpha}/2} \leq\frac{\sqrt{n}}{S}\left(\overline{X} - \beta \sqrt{\frac{\pi}{2}}\right) \leq {\zeta}_{1-{\alpha}/2}) \approx 1-{\alpha} \]
\[ P\!\left(-{\zeta}_{1-{\alpha}/2} \frac{s}{\sqrt{n}} \;\leq\; \overline{X} - \beta \sqrt{\tfrac{\pi}{2}} \;\leq\; {\zeta}_{1-{\alpha}/2} \frac{s}{\sqrt{n}}\right) \;\approx\; 1-\alpha. \]
\[ P\!\left(\overline{X} - {\zeta}_{1-{\alpha}/2}\frac{s}{\sqrt{n}} \;\leq\; \beta \sqrt{\tfrac{\pi}{2}} \;\leq\; \overline{X} + {\zeta}_{1-{\alpha}/2}\frac{s}{\sqrt{n}}\right) \;\approx\; 1-\alpha. \]
\[ P\!\left(\frac{\overline{X} - {\zeta}_{1-{\alpha}/2}\tfrac{s}{\sqrt{n}}}{\sqrt{\pi/2}} \;\leq\; \beta \;\leq\; \frac{\overline{X} + \zeta_{1-\alpha/2}\tfrac{s}{\sqrt{n}}}{\sqrt{\pi/2}}\right) \;\approx\; 1-\alpha. \]
Por lo tanto, obtenemos que \[ IC_{1-\alpha}(\beta) \;=\; \left[ \frac{\overline{X} - \zeta_{1-\alpha/2}\tfrac{s}{\sqrt{n}}}{\sqrt{\pi/2}}, \;\; \frac{\overline{X} + \zeta_{1-\alpha/2}\tfrac{s}{\sqrt{n}}}{\sqrt{\pi/2}} \right]. \]
De este modo, tomando \(\alpha = 0.02\), tenemos que \(\zeta_{1-\alpha/2} = \zeta_{0.99} \approx 2.326\), por lo que: \[ IC_{0.98}(\beta) \;=\; \left[ \frac{\overline{X} - 2.326 \cdot \tfrac{s}{\sqrt{n}}}{\sqrt{\pi/2}}, \;\; \frac{\overline{X} + 2.326 \cdot \tfrac{s}{\sqrt{n}}}{\sqrt{\pi/2}} \right]. \]
Ahora, probemos que \[ \frac{\sqrt{n}}{2\beta^2} \left( \overline{V} - 2\beta^2 \right) \approx N(0,1). \]
Sea \[ S_n = V_1 + V_2 + \dots + V_n. \] Sabemos que \[ \mathbb{E}(V_i) = \mathbb{E}(X_i^2) = \beta^2 \, 2^{2/2} \, \Gamma\!\left(1 + \frac{2}{2}\right) = 2\beta^2 \, \Gamma(2) = 2\beta^2 = \mu. \]
Además, \[ \mathbb{E}(V_i^2) = \beta^4 \, 2^{4/2} \, \Gamma\!\left(1 + \frac{4}{2}\right) = 4\beta^4 \, \Gamma(3) = 8\beta^4. \]
Así, \[ \sigma = \sqrt{\mathbb{E}(V_i^2) - \mathbb{E}(V_i)^2} = \sqrt{8\beta^4 - 4\beta^4} = \sqrt{4\beta^4} = 2\beta^2. \]
De este modo, \[ \frac{\sqrt{n}}{2\beta^2}\left(\overline{V} - 2\beta^2\right) = \frac{\sqrt{n}}{\sigma} \left(\overline{V} - \mu\right) = \frac{\sqrt{n}}{\sigma} \left(\frac{S_n - n\mu}{n}\right) = \frac{S_n - n\mu}{\sqrt{n}\,\sigma}. \]
Por el Teorema del Límite Central, \[ \frac{S_n - n\mu}{\sqrt{n}\,\sigma} \;\; \longrightarrow \;\; N(0,1), \quad \text{cuando } n \to \infty. \]
Luego, \[ \frac{S_n - n\mu}{\sqrt{n}\,\sigma} \approx N(0,1). \]
Ahora, deduzcamos un intervalo del \(1-\alpha\) de confianza para \(\beta\).
Del mismo modo que el intervalo calculado anteriormente, tenemos que
\[ \mathbb{P}\!\left( - \zeta_{1-\alpha/2} \leq \frac{\sqrt{n}}{2\beta^2}\left(\overline{V} - 2\beta^2\right) \leq \zeta_{1-\alpha/2} \right) \approx 1-\alpha. \]
\[ \mathbb{P}\!\left( - \frac{2\beta^2}{\sqrt{n}} \, \zeta_{1-\alpha/2} \;\;\leq\;\; \overline{V} - 2\beta^2 \;\;\leq\;\; \frac{2\beta^2}{\sqrt{n}} \, \zeta_{1-\alpha/2} \right) \approx 1-\alpha. \]
\[ \mathbb{P}\!\left( \overline{V} - \frac{2\beta^2}{\sqrt{n}} \, \zeta_{1-\alpha/2} \;\;\leq\;\; 2\beta^2 \;\;\leq\;\; \overline{V} + \frac{2\beta^2}{\sqrt{n}} \, \zeta_{1-\alpha/2} \right) \approx 1-\alpha. \]
\[ \mathbb{P}\!\left( \frac{\overline{V}}{2} - \frac{2\beta^2}{2\sqrt{n}} \, \zeta_{1-\alpha/2} \;\;\leq\;\; \beta^2 \;\;\leq\;\; \frac{\overline{V}}{2} + \frac{2\beta^2}{2\sqrt{n}} \, \zeta_{1-\alpha/2} \right) \approx 1-\alpha. \]
Finalmente, un intervalo de confianza para \(\beta\) es:
\[ IC_{1-\alpha}(\beta)= \left[ \sqrt{\tfrac{1}{2}\left(\overline{V} - \tfrac{2\beta^2}{\sqrt{n}} \, \zeta_{1-\alpha/2}\right)}, \;\; \sqrt{\tfrac{1}{2}\left(\overline{V} + \tfrac{2\beta^2}{\sqrt{n}} \, \zeta_{1-\alpha/2}\right)} \right]. \]
Sea \(X_1,\dots,X_n\) una muestra
i.i.d. de una variable con f.d.a \(F\)
de la Rayleigh(\(\beta\)) .
Sean
\[ U_i := F(X_i), \] Entonces \(U_1,\dots,U_n\) son i.i.d. \(\operatorname{Unif}(0,1)\).
Además, tenemos que
\[ F\!\big(X_{(k)}\big) = U_{(k)}. \]
Sea \(U_{(k)}\) la \(k\)-ésima variable ordenada de una muestra
i.i.d. \(\operatorname{Unif}(0,1)\).
La función de distribución acumulada esta dada por \[
\begin{aligned}
\Pr\big(U_{(\frac{n}{2})} \le t\big)
&= \Pr\big(\text{al menos $\frac{n}{2}$ de los $U_i$ son } \le
t\big) \\
&= \sum_{j=\frac{n}{2}}^{n} \binom{n}{j} t^{j}(1-t)^{n-j}.
\end{aligned}
\] Al derivar para obtener la densidad, el resultado es \[
f_{U_{(\frac{n}{2})}}(t) = \frac{n!}{(\frac{n}{2}-1)!(n-\frac{n}{2})!}\,
t^{\frac{n}{2}-1}(1-t)^{n-\frac{n}{2}}.
\] Por lo que obtenemos que
\[ F\big(X_{(n/2)}\big) \sim \mathrm{Beta}\!\left(\tfrac{n}{2},\, \tfrac{n}{2}+1\right). \]
Sean \(c_\ell\) y \(c_u\) los cuantiles \(\alpha/2\) y \(1-\alpha/2\), respectivamente, de la distribución Beta obtenida; entonces, tenemos que
\[ \Pr\!\left(c_\ell \leq F(X_{(n/2)}) \leq c_u \right) = 1-\alpha. \]
\[ \Pr\!\left(F^{-1}(c_\ell) \leq X_{(n/2)} \leq F^{-1}(c_u)\right) = 1-\alpha. \]
Sin embargo, como calculamos en la primera parte,
\[ F^{-1}(x) = \beta \sqrt{-2\ln(1-x)} \]
así,
\[ \beta \sqrt{-2\ln(1-c_\ell)} \;\leq\; X_{(n/2)} \;\leq\; \beta \sqrt{-2\ln(1-c_u)}. \]
Por lo tanto, \[ \Pr\!\left(\, \frac{X_{(n/2)}}{\sqrt{-2\ln(1-c_u)}} \;\leq\; \beta \;\leq\; \frac{X_{(n/2)}}{\sqrt{-2\ln(1-c_\ell)}} \right) = 1-\alpha. \]
Finalmente, un intervalo del \(1-\alpha\) de confianza para \(\beta\) es \[ IC_{1-\alpha}(\beta) \;=\; \left[ \frac{X_{(n/2)}}{\sqrt{-2\ln(1-c_u)}},\; \frac{X_{(n/2)}}{\sqrt{-2\ln(1-c_\ell)}} \right]. \]
Encontremos ahora otro procedimiento de intervalo de confianza para β.
Sea \(V_i = X_i^4\) y \(\overline{V} = \frac{1}{n}\sum_{i=1}^n V_i\). Sabemos que \[ \mathbb{E}[V_i]=\mathbb{E}[X^4]=8\beta^4, \qquad \mathbb{E}[V_i^2]=\mathbb{E}[X^8]=384\,\beta^8; \] así, \[ \operatorname{Var}(V_i)=384\beta^8-(8\beta^4)^2=320\beta^8. \]
Como resultafo, \(\sigma=\sqrt{320}\,\beta^4=8\sqrt{5}\,\beta^4\).
Por el Teorema del Límite Central, \[ \frac{\sqrt{n}\,(\overline{V}-8\beta^4)}{8\sqrt{5}\,\beta^4}\xrightarrow{d} N(0,1). \]
Al igual que en los anteriores casos, para el cuantil de \(1-\alpha\)
\[ \Pr\!\Big(-\zeta_{1-\alpha/2}\le \frac{\sqrt{n}(\overline{V}-8\beta^4)}{8\sqrt{5}\,\beta^4}\le \zeta_{1-\alpha/2}\Big)\approx 1-\alpha. \]
Obtenemos que
\[ \Pr\!\left(\frac{\overline{V}}{8 + \dfrac{8\zeta_{1-\alpha/2}\sqrt{5}}{\sqrt{n}}}\le \beta^4 \le \frac{\overline{V}}{8 - \dfrac{8\zeta_{1-\alpha/2}\sqrt{5}}{\sqrt{n}}}\right)\approx 1-\alpha. \] Finalmente, llegamos a que el intervalo es
\[ IC_{1-\alpha}(\beta)\approx\left[ \left(\frac{\overline{V}}{8 + \tfrac{8\zeta_{1-\alpha/2}\sqrt{5}}{\sqrt{n}}}\right)^{1/4},\; \left(\frac{\overline{V}}{8 - \tfrac{8\zeta_{1-\alpha/2}\sqrt{5}}{\sqrt{n}}}\right)^{1/4} \right]. \]
# Generamos la muestra
set.seed(123)
B <- 50000
beta <- 2.5
n_vals <- c(50, 100, 200, 500)
alpha <- 0.02
# Función generadora de Rayleigh(β) que calculamos en la primera parte
rRayleigh <- function(n, beta){
sqrt(-2*beta^2 * log(runif(n)))
}
# Método 1
IC1_mean <- function(x, beta, alpha){
n <- length(x)
s <- sd(x)
xbar <- mean(x)
mu <- beta*sqrt(pi/2)
zeta <- qnorm(1 - alpha/2)
L <- (xbar - zeta*s/sqrt(n)) / sqrt(pi/2)
U <- (xbar + zeta*s/sqrt(n)) / sqrt(pi/2)
c(max(L,0), U) # límite inferior no negativo
}
# Método 2
IC2_square <- function(x, beta, alpha){
n <- length(x)
V <- x^2
Vbar <- mean(V)
zeta <- qnorm(1 - alpha/2)
L <- sqrt((Vbar - zeta*2*beta^2/sqrt(n)) / 2)
U <- sqrt((Vbar + zeta*2*beta^2/sqrt(n)) / 2)
c(max(L,0), U)
}
# Método 3
IC3_median <- function(x, beta, alpha){
n <- length(x)
m <- sort(x)[n/2] # mediana muestral (n par)
cl <- qbeta(alpha/2, n/2, n/2+1)
cu <- qbeta(1-alpha/2, n/2, n/2+1)
# Invertir la CDF de Rayleigh
L <- m / sqrt(-2*beta^2*log(1-cu))
U <- m / sqrt(-2*beta^2*log(1-cl))
c(max(L,0), U)
}
# Método 4
IC4_fourth <- function(x, beta, alpha){
n <- length(x)
V <- x^4
Vbar <- mean(V)
mu <- 8*beta^4
sigma <- sqrt(384*beta^8 - mu^2) # Var(X^4)
zeta <- qnorm(1 - alpha/2)
L <- ((Vbar - zeta*sigma/sqrt(n))/8)^(1/4)
U <- ((Vbar + zeta*sigma/sqrt(n))/8)^(1/4)
c(max(L,0), U)
}
# Monte Carlo
resultados <- list()
# Calculamos la cobertura y longitud
for (n in n_vals){
cobertura <- matrix(0, nrow=B, ncol=4)
longitudes <- matrix(0, nrow=B, ncol=4)
for (b in 1:B){
x <- rRayleigh(n, beta)
ICs <- list(
IC1_mean(x, beta, alpha),
IC2_square(x, beta, alpha),
IC3_median(x, beta, alpha),
IC4_fourth(x, beta, alpha)
)
for (m in 1:4){
L <- ICs[[m]][1]
U <- ICs[[m]][2]
longitudes[b,m] <- U - L
cobertura[b,m] <- (beta >= L & beta <= U)
}
}
resultados[[as.character(n)]] <- data.frame(
Metodo = paste0("IC", 1:4),
Cobertura = colMeans(cobertura),
SE_cobertura = sqrt(colMeans(cobertura)*(1-colMeans(cobertura))/B),
Longitud = colMeans(longitudes),
SE_longitud = apply(longitudes, 2, sd)/sqrt(B)
)
}
resultados
## $`50`
## Metodo Cobertura SE_cobertura Longitud SE_longitud
## 1 IC1 0.97418 0.0007092718 0.8548628 0.0004080274
## 2 IC2 0.97946 0.0006343202 0.8414549 0.0002873662
## 3 IC3 0.00000 0.0000000000 0.4937307 0.0002248378
## 4 IC4 NA NA NaN NA
##
## $`100`
## Metodo Cobertura SE_cobertura Longitud SE_longitud
## 1 IC1 0.97826 0.0006521867 0.6063324 0.0002042205
## 2 IC2 0.98016 0.0006236405 0.5880247 0.0001359414
## 3 IC3 0.00000 0.0000000000 0.3421741 0.0001106240
## 4 IC4 NA NA NaN NA
##
## $`200`
## Metodo Cobertura SE_cobertura Longitud SE_longitud
## 1 IC1 0.97878 0.0006445108 0.4292165 1.017698e-04
## 2 IC2 0.97970 0.0006306807 0.4134845 6.659146e-05
## 3 IC3 0.00000 0.0000000000 0.2396581 5.481042e-05
## 4 IC4 0.97900 0.0006412332 0.4845969 2.847859e-04
##
## $`500`
## Metodo Cobertura SE_cobertura Longitud SE_longitud
## 1 IC1 0.97906 0.0006403361 0.2717784 4.058870e-05
## 2 IC2 0.97986 0.0006282417 0.2606426 2.607429e-05
## 3 IC3 0.00000 0.0000000000 0.1506783 2.176211e-05
## 4 IC4 0.98066 0.0006158890 0.2964341 1.021989e-04
De esta tabla, comparamos los métodos IC1, IC2, IC3, IC4. En este contexto, la cobertura hace referencia a la probabilidad estimada de contener el parámetro; la longitud determina el ancho del intervalo y los SE_cobertura y SE_longitud son los respectivos errores estándar (incertidumbre en los estimadores producidos).
De los resultados obtenidos, podemos evidenciar que los mejores métodos fueron IC1 e IC2, ya que para cada valor de n la probabilidad estimada de contener el parámetro fue bastante alta y la incertidumbre fue bastante pequeña; en adición, la longitud de los intervalos era mínima e incluso se redujo bastante cuando n tomaba valores más grandes. Por parte del método IC3, se evidenció que no fue un método útil, debido a que para cada valor de n la probabilidad de estimada de contener el parámetro fue cero, además que la longitud de los intervalos siempre fue mayor a cero. Por parte del método IC4, evidenciamos que para valores de n pequeños no fue un método práctico, ya que con este método no se lograron calcular valores para la cobertura ni longitud; sin embargo, para valores más grandes de n, el método funcionó bastante bien generando una probailidad de contener el parámetro muy alta con logitudes pequeñas y incertidumbres bastante bajas para ambos parámetros. En general, para cada método que se comportó de buena manera, se evidencio que a medida que n tomaba valores más grandes la cobertura era mayor, la longitud menor y las respectivas incertidumbre eran más bajas.
Con respecto a la complejidad computacional, los métodos IC1 e IC2 no generaron tanto costo para ser calculado. Por otra parte, el método IC4 requería un poco más de costo en comparación con los anteriores. Finalmente, el método IC3 requería bastante costo, esto resultó en que los valores obtenidos en la tabla no fueran los valores esperados. La variación en el peso de estos cálculos puede deberse a la complejidad en las ecuaciones para ser determinados, ya sea por los exponentes, raices o potencias que pudieran contener.
A partir del desarrollo del proyecto y de los resultados obtenidos, se pueden resaltar los siguientes puntos importantes:
El algoritmo 1 para obtener una muestra con la distribución de Rayleigh a partir de una muestra uniforme, resultó ser bastante efectivo. A través del q-q plot logramos confirmar esta afirmación commparándola con la distribución teórica.
Los diferentes métodos de obtención de intervalos de confianza fueron, en general, útiles para la estimación del parámetro \(\beta\). Aún así, se evidenció una mejor estimación para tamaños muestrales mayores.
A tráves de las simulaciones de Monte Carlo, logramos comparar los diferentes métodos de estimación. Además, observamos que a medida que el tamaño muestral aumenta, la cobertura, longitud e incertidumbre se comportan de mejor manera
Por medio del proyecto, logramos profundizar en las propiedades de la distribución de Rayleigh y hacer comparaciones entre distribuciones muestrales y teóricas, además de los diferentes métodos y efectividad de cada uno de ellos.