Ejercicio 1, Práctico 2

Parte a

Calculemos el valor esperado de \(X \sim Γ(α, β)\) y su varianza. El valor esperado de una v.a. con pdf \(f(x)\) se calcula como la integral

\[ E(X) = \int_{-\infty}^\infty x f(x) dx, \] y dado que \(f(x) = f(x \vert \alpha, \beta) = (\Gamma(\alpha)\beta^\alpha)^{-1} x^{\alpha-1}e^{-x/\beta}\), para \(0 < x < \infty\) y \(\alpha, \beta > 0\), tenemos que

\[ E(X) = \dfrac{1}{\Gamma(\alpha)\beta^\alpha} \int_{0}^\infty x^\alpha e^{-x/\beta}dx. \] Con el c.v. \(u = x/β\), la integral queda:

\[ \int_{0}^\infty x^\alpha e^{-x/\beta}dx = β^{α+1}\int_{0}^\infty u^α e^{-u}du = β^{α+1} Γ(α+1). \]

Sustituyendo en \(E(x)\), y recordando que \(Γ(z+1) = zΓ(z)\) para todo \(z\) con \(ℜ(z) > -1\), concluimos que:

\[ E(x) = \dfrac{1}{Γ(α)β^α}β^{α+1} Γ(α+1) = αβ. \]

Calculemos la varianza de \(X\). Con \(E(X^2)\) nos encontramos con la integral

\[ \int_{0}^\infty x^{\alpha+1} e^{-x/\beta}dx = β^{α+2}\int_{0}^\infty u^{α+1} e^{-u}du = β^{α+2} Γ(α+2). \]

Por tanto, \(E(X^2) = (β^{α+2} Γ(α+2))/(β^{α} Γ(α)) = β^2(α+1)α\), y

\[ Var(X) = E(X^2) - E(X)^2 = (α+1)αβ^2 - α^2β^2 = αβ^2. \]

Parte b

Recordemos la noción de función de verosimilitud.

Def. (función de verosimilitud). Sea \(f(x\vert θ)\) la función pdf conjunta o mdf de la muestra \(X = (X_1, \ldots, X_n)\). Entonces, dado que \(X = x\) fue observado, la función de \(\theta\) definida como

\[ L(\theta \vert x) = f(x \vert \theta) \] se llama función de verosimilitud.

La conexión entre la estimación puntual de un parámetro y la función de máxima verosimilitud, \(L(\theta \vert x)\), está comprendida por la siguiente definición de estimador máximo verosímil (MLE).

Def. (MLE). Para cada dato muestral \(x\), sea \(\hat{\theta}\) el valor de un parámetro en el cual \(L(\theta \vert x)\) alcanza su máximo como función de \(\theta\), dejando a \(x\) fijo. Un estimador de máxima verosimilitud (MLE) del parámetro \(\theta\) basado en la muestra \(X\) es \(\hat{\theta}(X)\).

En esta parte asumimos que el parámetro \(α\) es conocido y deseamos realizar la estimación del parámetro \(β\) utilizando el MLE asociado. La función de verosimilitud es

\[ L(\beta\vert x) = ∏_{i=1}^n \dfrac{1}{Γ(α)β^{α}}x_i^{α-1}e^{-x_i/β} = \dfrac{1}{Γ(α)^n β^{nα}} \left[ ∏_{i=1}^n x_i \right]^{α-1} e^{-∑_i x_i/β}. \]

Es conveniente trabajar con la log-verosimilitud. Como el logaritmo es una función monótona, a los efectos de hallar el máximo de \(L\) podemos calcularlo utilizando la log-verosimilitud.

\[ \log L(β \vert x) = -n\log \Gamma(\alpha) - n\alpha \log \beta + (\alpha-1)\sum_i \log x_i - \beta^{-1}\sum_i x_i. \] Por lo tanto, \[ \dfrac{\partial \log L}{\partial \beta} = -\dfrac{n \alpha}{\beta} + \beta^{-2}\sum_i x_i. \] Haciendo esta derivada parcial igual a cero, y resolviendo para \(\beta\), encontramos que \[ \hat{\beta} = \dfrac{\sum_i x_i}{n\alpha} = \dfrac{\bar{X}}{\alpha}. \]

Verifiquemos que efectivamente se trata de un máximo:

\[ \dfrac{\partial^2 \log L}{\partial \beta^2}_{\vert \beta = \hat{\beta}} = \left[ \frac{n\alpha}{\beta^2} - \frac{2\sum_i x_i}{\beta^3}\right]_{\vert \beta = \hat{\beta}} = -\dfrac{(n\alpha)^3}{(\sum_i x_i)^2} < 0. \]

Dado que \(\hat{\beta}\) es el único punto donde la derivada de \(L\) es cero y es un máximo local, entonces es un máximo global. Deducimos que \(\hat{\beta}\) es el MLE.

Parte c

En esta parte, la función de verosimilitud es

\[ L(α, β \vert x) = \dfrac{1}{Γ(α)^n β^{nα}} \left[ ∏_{i=1}^n x_i \right]^{α-1} e^{-∑_i x_i/β}. \]

Es decir la misma que en la parte (a), salvo que ahora tanto \(\beta\) como \(\alpha\) se consideran variables.

El cálculo analítico de \(\hat{\alpha}\) y de \(\hat{\beta}\) no se puede hacer utilizando funciones elementales, ya que intervienen derivadas de la función gamma. Estas derivadas se describen con la función especial polygamma. Es decir, los valores analíticos de \(\hat{\alpha}\) y de \(\hat{\beta}\) nos quedarían expresados en términos de ciertas funciones especiales.

El enfoque alternativo sugerido en este ejercicio es realizar el cálculo de los MLE numéricamente, con la salvedad siguiente. Conocemos el máximo de \(L\) para cada valor fijo de \(α\), e inyectamos el maximizador (\(β = β(α)\)) en la fórmula de \(L(α, β)\). De esta forma hemos transformado un problema en principio bidimensional a un problema de optimización esencialmente unidimensional.

Entonces por la parte (b), tenemos una fórmula para \(\hat{\beta}\) como función de \(\hat{\alpha}\). Por tanto, podemos sustituir dicha fórmula en \(L\):

\[ L(α, β=\bar{X}/\alpha \vert x) = \dfrac{1}{Γ(α)^n (\bar{X}/\alpha)^{nα}} \left[ ∏_{i=1}^n x_i \right]^{α-1} e^{-n\alpha}. \]

A continuación escribimos un programa en R para calcular el máximo de esta función dados los datos del problema.

El comando c de R se utiliza para concatenar argumentos:

# ingresar datos muestrales
x = c(22.0, 23.9, 20.9, 23.8, 25.0, 24.0, 21.7, 23.8, 22.8, 23.1, 23.1, 23.5, 23.0, 23.0)

Chequeo el primer elemento y la cantidad de datos en el vector x:

> x[1]
[1] 22

> n = length(x); n
[1] 14

Como x tiene números relativamente "grandes", conviene trabajar con la log-verosimilitud. Nótese que, por ejemplo

> gamma(20)
[1] 1.216451e+17

> log(gamma(20))
[1] 39.33988

> lgamma(20)
[1] 39.33988

Aquí lgamma es la función especial que corresponde al logaritmo natural del valor absoluto de la función gamma. Conviene trabajar con esta función ya que es más precisa que primero calcular gamma y luego tomar el logaritmo (ver ejemplos más abajo).

La log-verosimilitud es

\[ \log L(α, β=\bar{X}/α \vert x) = -n\log \Gamma(\alpha) - n\alpha \log (\bar{X}/\alpha) + (\alpha-1)\sum_i \log x_i - n\alpha. \]

La media muestral se calcula con la función mean:

> xmean = mean(x)
> xmean
[1] 23.11429

Primero ingresamos la función \(\log L(α, β=\bar{X}/α \vert x)\):

logL <- function(alpha){
  -n * lgamma(alpha) - n * alpha * log(xmean/alpha) + (alpha-1) * sum(log(x)) - n * alpha
}

Evaluamos algunos valores:

> logL(10)
[1] -40.95251
> logL(100)
[1] -25.95473
> logL(500)
[1] -20.12502

A modo de comparación, si hubiéramos utilizado directamente la función \(L\), no podríamos obtener estos valores, ya que escapamos del rango de valores admisibles de las funciones especiales que se utilizan en su definición:

> L <- function(alpha){
  prod(x)^(alpha-1) * exp(-n*alpha) / gamma(alpha)^n / (sum(x)/(n*alpha))^(alpha * n)
}

> L(10)
[1] 1.638892e-18
> L(100)
[1] NaN
> L(500)
[1] NaN
Warning message:
In L(500) : value out of range in 'gammafn'

Este ejemplo demuestra la importancia de considerar el dominio de validez de los métodos numéricos que utilizamos, para obtener valores correctos.

Continuando con la función log-verosimilitud, queremos hallar el máximo respecto a \(α\). Utilizaremos la función de R optimize para realizar dicho cálculo. En la documentación vemos que por defecto esta función hace una minimización. Una alternativa es trabajar con la minimización de \(-\log L\).

Por defecto, optimize utiliza un metodo de "line search" que necesita dos puntos para acotar el intervalo, tomemos 1 y 1000 con una tolerancia de una parte en un millon como criterio de terminacion:

> minus_logL <- function(alpha){
  n * lgamma(alpha) + n * alpha * log(xmean/alpha) - (alpha-1) * sum(log(x)) + n * alpha
}

> optimize(minus_logL, c(1, 1e3), tol=1e-6)
$minimum
[1] 514.3353

$objective
[1] 20.12224

> mean(x) / 514.3353
[1] 0.04494011

Concluimos que \(\hat{\alpha} = 514\) y \(\hat{\beta} = 0.0449\).

Nota. Inspeccionando la documentación de optimize vemos que tiene un argumento maximum que es por defecto falso (FALSE):

> optimize
function (f, interval, ..., lower = min(interval), upper = max(interval), 
    maximum = FALSE, tol = .Machine$double.eps^0.25) 
{
    if (maximum) {
        val <- .External2(C_do_fmin, function(arg) -f(arg, ...), 
            lower, upper, tol)
        list(maximum = val, objective = f(val, ...))
    }
    else {
        val <- .External2(C_do_fmin, function(arg) f(arg, ...), 
            lower, upper, tol)
        list(minimum = val, objective = f(val, ...))
    }
}
<bytecode: 0x1123912a8>
<environment: namespace:stats>

En nuestro caso, como interesa el máximo, podemos trabajar directamente con logL utilizando el comando optimize(logL, c(1, 1e3), maximum=TRUE, tol=1e-6).