markdown

Introducción

El objetivo de este proyecto es analizar diferentes aspectos de la convergencia del Estimador de Máxima Verosimilitud (EMV) biparamétrico del parámetro \(\theta = (\alpha, \beta)\) en el modelo Gamma, así como calcular los porcentajes de pertenecia al conjunto de confianza elíptico que se obtiene a partir del Teorema del Límite Central (TLC) aplicado al EMV.

Para ello, vamos a implementar el método de Newton–Raphson para estimar los parámetros \((\alpha, \beta)\) a partir de datos simulados bajo la distribución Gamma. Los valores iniciales del proceso serán determinados mediante los estimadores de momentos, lo cual garantiza una buena proximidad al EMV y acelera la convergencia del algoritmo.

A partir de esto, vamos a calcular:

La velocidad de convergencia del método de Newton–Raphson para distintos tamaños muestrales \(n\).
El comportamiento empírico de los conjuntos de confianza elípticos construidos con base en la matriz de información de Fisher.
Los porcentajes de pertenencia del parámetro verdadero \(\theta_0\) en los conjuntos de confianza y las áreas promedio de las elipses asociadas
La verificación de la convergencia

\[ n \cdot \mathrm{Cov}(\hat{\theta}) \approx I_F^{-1}(\theta_0). \]

Todo esto teniendo en cuenta los diferentes tamaños muestrales \(n\).

Nota: En el código de python los datos e información varian cada vez que se ejecuta el código. Los datos en este proyecto están determinandos por una ejecución especifica del código.

Estimación por Newton-Raphson

Calculemos \(l(\theta)\), con \(\theta = (\alpha, \beta)\), entonces

\[ L(\theta) = \prod_{i=1}^{n} \frac{1}{\Gamma(\alpha)\beta^{\alpha}} x_i^{\alpha - 1} e^{-x_i / \beta} = \left( \frac{1}{\Gamma(\alpha)\beta^{\alpha}} \right)^n \prod_{i=1}^{n} x_i^{\alpha - 1} e^{-x_i / \beta}. \]

De este modo,

\[ l(\theta) = \ln(L(\theta)) = n \ln \left( \frac{1}{\Gamma(\alpha)\beta^{\alpha}} \right) + (\alpha - 1) \sum_{i=1}^{n} \ln(x_i) - \frac{ln(e)}{\beta} \sum_{i=1}^{n} x_i, \] \[ = -n \left( \ln(\Gamma(\alpha)) + \alpha \ln(\beta) \right) + (\alpha - 1) \sum_{i=1}^{n} \ln(x_i) - \frac{1}{\beta} \sum_{i=1}^{n} x_i. \]

Calculando \(\dfrac{\partial l}{\partial \alpha}\) y \(\dfrac{\partial l}{\partial \beta}\), obtenemos

\[ \frac{\partial l}{\partial \alpha} = -n \left( \psi(\alpha) + \ln(\beta) \right) + \sum_{i=1}^{n} \ln(x_i), \]

\[ \frac{\partial l}{\partial \beta} = -\frac{n \alpha}{\beta} + \frac{1}{\beta^2} \sum_{i=1}^{n} x_i, \] donde \(\psi(\alpha) = \frac{\Gamma'(\alpha)}{\Gamma(\alpha)}.\)

Dado \(x_0 = (\alpha_0, \beta_0)^{T}\) vector inicial, entonces utilizando el método de Newton-Raphson obtenemos la sucesión

\[ x_k = x_{k-1} - [DF(x_{k-1})]^{-1} f(x_{k-1}), \] tal que \[ f(x_k) = f(\alpha_k, \beta_k) = \left( \frac{\partial l}{\partial \alpha}(\alpha_k, \beta_k), \, \frac{\partial l}{\partial \beta}(\alpha_k, \beta_k) \right). \]

Luego,

\[ DF = H = \begin{pmatrix} -n \psi'(\alpha) & -\dfrac{n}{\beta} \\ -\dfrac{n}{\beta} & \dfrac{n \alpha}{\beta^2} -\dfrac{2}{\beta^3} \sum_{i=1}^{n} x_i \end{pmatrix}. \]

Vamos a generar el código en python, por lo que para implementar el método de Newton-Raphson, no vamos a calcular la matriz inversa de \(H\); sino que vamos a utilizar la función “linalg.solve” de la libreria “Numpy” para calcular el \([DF(x_k)]^{-1}f(x_k)\).

Por otra parte, para estimar los valores iniciales \((\alpha_0, \beta_0)\) para el método de Newton-Raphson, tenemos que

\[ \mu = \alpha \beta \quad \text{y} \quad \sigma^2 = \alpha \beta^2. \]

Luego,

\[ \mu \beta = \alpha \beta^2 = \sigma^2, \quad \text{es decir,} \quad \beta = \frac{\sigma^2}{\mu}. \]

Además,

\[ \alpha = \frac{\mu}{\beta} = \frac{\mu^2}{\sigma^2}. \]

Como los estimadores de \(\mu\) y \(\sigma^2\) son \(\bar{X}\) y \(s^2\), respectivamente, entonces

\[ \tilde{\alpha} = \frac{\bar{X}^2}{s^2} \quad \text{y} \quad \tilde{\beta} = \frac{s^2}{\bar{X}}. \]

El link de Google Colab donde ejecutamos el código con las muestras, método de N-R y gráficas es el siguiente: https://colab.research.google.com/drive/1H25sHU00KPQZr5nEjkYDYnfK-kbo6KRA?usp=sharing

Gráficas

Las gráficas y boxplots obtenidos son los siguientes:

Además, como se puede visualizar en el código, la cantidad de iteraciones obtenidas son:

“el número aproximado de iteraciones que requiere N-R para converger con n=100 es 2.”

“el número aproximado de iteraciones que requiere N-R para converger con n=200 es 3.”

“el número aproximado de iteraciones que requiere N-R para converger con n=500 es 3.”

“el número aproximado de iteraciones que requiere N-R para converger con n=1000 es 2.”

“el número aproximado de iteraciones que requiere N-R para converger con n=2000 es 2.”

Esto nos dice que la cantidad de iteraciones que necesita el proceso de Newton-Raphson es bastante pequeña, en general menor a 4. Esto se debe a que los valores iniciales de este proceso están bastante cerca del EMV.

Por otra parte, de las gráficas podemos evidenciar que, para cada tamaño muestral, las estimaciones se concentran cerca del punto \((\alpha=2.5 ,\beta= 4.6)\) tanto en las graficas \(\alpha \times \beta\) como en los boxplots. Sin embargo, mientras el tamaño de la muestra es más grande, las estimaciones son menos dispersas; es decir, todas las estimaciones pertenecen a un rango más cercano del EMV.

Conjuntos de confianza

Primero, encontremos la matriz \(I_F(\theta)\).

Sea \[ l_1 = \ln \left( \frac{1}{\Gamma(\alpha)\beta^\alpha} x^{\alpha-1} e^{-x/\beta} \right) = -\ln(\Gamma(\alpha)) - \alpha \ln(\beta) + (\alpha - 1)\ln(x) - \frac{x}{\beta}. \]

Entonces,

\[ I_F(\theta) = -\mathbb{E} \begin{pmatrix} \frac{\partial^2 l_1}{\partial \alpha^2} & \frac{\partial^2 l_1}{\partial \alpha \partial \beta} \\[6pt] \frac{\partial^2 l_1}{\partial \beta \partial \alpha} & \frac{\partial^2 l_1}{\partial \beta^2} \end{pmatrix}_{(\alpha,\beta)} \]

\[ = -\mathbb{E} \begin{pmatrix} \frac{\partial}{\partial \alpha} \left[-\psi(\alpha) - \ln(\beta) + \ln(x)\right] & \frac{\partial}{\partial \alpha} \left[-\frac{\alpha}{\beta} + \frac{x}{\beta^2}\right] \\[6pt] \frac{\partial}{\partial \beta} \left[-\psi(\alpha) - \ln(\beta) + \ln(x)\right] & \frac{\partial}{\partial \beta} \left[-\frac{\alpha}{\beta} + \frac{x}{\beta^2}\right] \end{pmatrix}_{(\alpha,\beta)}. \]

\[ = -\mathbb{E} \begin{pmatrix} -h(\alpha) & -\frac{1}{\beta} \\[6pt] -\frac{1}{\beta} & \frac{\alpha}{\beta^2} - \frac{2x}{\beta^3} \end{pmatrix}_{(\alpha,\beta)}. \]

Dado que \(\mathbb{E}(X) = \alpha \beta\), entonces

\[ I_F(\theta) = -\begin{pmatrix} -h(\alpha) & -\frac{1}{\beta} \\[6pt] -\frac{1}{\beta} & \frac{\alpha}{\beta^2}- \frac{2 \alpha \beta}{\beta^3} \end{pmatrix} = \begin{pmatrix} h(\alpha) & \frac{1}{\beta} \\[6pt] \frac{1}{\beta} & \frac{\alpha}{\beta^2} \end{pmatrix}. \]

Por el Teorema del Límite Central (TLC), para \(\hat{\theta}\) tenemos que

\[ \sqrt{n} \, I_F(\hat{\theta})^{1/2} (\hat{\theta} - \theta_0) \xrightarrow{d} N(0, I), \]

donde \(I\) es la matriz identidad.

De esto, tenemos que

\[ n (\hat{\theta} - \theta_0)^T I_F(\hat{\theta}) (\hat{\theta} - \theta_0) \xrightarrow{d} \chi^2_2. \]

De este modo, sea \(\rho\) el cuantil del 95% de \(\chi^2_2\), entonces el conjunto

\[ C = \left\{ \theta : n(\hat{\theta} - \theta)^T I_F(\hat{\theta}) (\hat{\theta} - \theta) \leq \rho \right\} \]

tiene probabilidad del 95% de contener a \(\theta_0\).

De este modo, calculando el porcentaje de veces que \(\theta_0\) pertenece a \(C\) con los valores obtenidos en la primera parte, dados por el código en el link

https://colab.research.google.com/drive/1H25sHU00KPQZr5nEjkYDYnfK-kbo6KRA?usp=sharing

obtenemos,

n=100: porcentaje empírico = 91.20%

n=200: porcentaje empírico = 92.90%

n=500: porcentaje empírico = 94.70%

n=1000: porcentaje empírico = 94.10%

n=2000: porcentaje empírico = 95.30%

De los porcentajes obtenidos, podemos evidenciar que todos estan bastante cerca del 95%.Sin embargo, a medida que el tamaño muestral es más grande, podemos ver que el porcentaje se acerca aún más al 95%. De manera aproximada, a partir de n=500 se cumple de manera muy acertada la probabilidad de contención del valor \(\theta_0\)

Sea \(U\) una matriz \(2 \times 2\) ortogonal tal que \[ U^{\top} I_F(\hat{\theta}) U = D, \] donde \(D\) es la matriz diagonal con los valores propios de \(I_F(\hat{\theta})\).

Entonces, la condición de que \(\theta\) pertenezca a \(C\) es equivalente a \[ n (\hat{\theta} - \theta)^{\top} U D U^{\top} (\hat{\theta} - \theta) \leq \rho. \]

Esto nos dice que \(C\) es una elipse en la variable

\[ Y = U^{\top} (\hat{\theta} - \theta). \]

Calculando el área de esta elipse para cada muestra en el código que se encuentra en el link https://colab.research.google.com/drive/1H25sHU00KPQZr5nEjkYDYnfK-kbo6KRA?usp=sharing

tenemos que,

n=100: área promedio = 1.801

n=200: área promedio = 0.9096

n=500: área promedio = 0.3627

n=1000: área promedio = 0.182

n=2000: área promedio = 0.091

Convergencia de la covarianza

Tenemos que

\[ I_F^{-1}(\hat{\theta}) = \frac{1}{\det(I_F(\hat{\theta}))} \begin{pmatrix} \dfrac{\alpha}{\beta^2} & -\dfrac{1}{\beta} \\[6pt] -\dfrac{1}{\beta} & h(\alpha) \end{pmatrix} \]

\[ I_F^{-1}(\hat{\theta}) = \frac{1}{\dfrac{\alpha h(\alpha)}{\beta^2} - \dfrac{1}{\beta^2}} \begin{pmatrix} \dfrac{\alpha}{\beta^2} & -\dfrac{1}{\beta} \\[6pt] -\dfrac{1}{\beta} & h(\alpha) \end{pmatrix} \]

\[ I_F^{-1}(\hat{\theta}) = \frac{\beta^2}{\alpha h(\alpha) - 1} \begin{pmatrix} \dfrac{\alpha}{\beta^2} & -\dfrac{1}{\beta} \\[6pt] -\dfrac{1}{\beta} & h(\alpha) \end{pmatrix} \]

\[ I_F^{-1}(\hat{\theta}) = \begin{pmatrix} \dfrac{\alpha}{\alpha h(\alpha) - 1} & \dfrac{\beta}{1 - \alpha h(\alpha)} \\[8pt] \dfrac{\beta}{1 - \alpha h(\alpha)} & \dfrac{\beta^2 h(\alpha)}{\alpha h(\alpha) - 1} \end{pmatrix} \]

Por el Teorema del Límite Central (TLC) para \(\hat{\theta}\), la covarianza muestral debería aproximarse a la matriz

\[ {I_F^{-1}(\hat{\theta})}/{n} = \begin{pmatrix} \dfrac{\alpha}{n(\alpha h(\alpha) - 1)} & \dfrac{\beta}{n(1 - \alpha h(\alpha))} \\[8pt] \dfrac{\beta}{n(1 - \alpha h(\alpha))} & \dfrac{\beta^2 h(\alpha)}{n(\alpha h(\alpha) - 1)} \end{pmatrix} \]

Realicemos la comparación de la covarianza y \(I_F^{-1}(\hat{\theta})\) con los datos ya obtenidos en el código, calculando \(n \, \text{Cov}(\hat{\theta}_1, \ldots, \hat{\theta}_m)\) y la diferencia dada por la norma de Frobenius, \[ \| n \, \text{Cov}(\hat{\theta}_1, \ldots, \hat{\theta}_m) - I_F^{-1}(\hat{\theta}) \|. \]

A partir del código en el link

https://colab.research.google.com/drive/1H25sHU00KPQZr5nEjkYDYnfK-kbo6KRA?usp=sharing

obtenemos los resultados:

\(n = 100\): \[ n \cdot \text{Cov}(\hat{\theta}_1, ... , \hat{\theta}_m) = \begin{pmatrix} 11.90589971 & -20.87728174 \\ -20.87728174 & 46.34487264 \end{pmatrix} \]