El objetivo de esta segunda práctica es aprender a calcular los estimadores de máxima verosimilitud de los parámetros de un modelo paramétrico en la práctica, aunque no sea posible obtener una expresión matemática explícita de los mismos. El modelo paramétrico empleado es una gamma.
La distribución gamma de parámetros \(\alpha>0\) y \(\beta>0\), \(\Gamma(\alpha,\beta)\), tiene función de densidad \[ f(x;\alpha,\beta)= \frac{1}{\Gamma(\alpha)\beta^\alpha} \, x^{\alpha-1} e^{-x/\beta} \; , \quad 0<x<\infty \, , \] donde \[ \Gamma(\alpha) = \int_0^\infty t^{\alpha-1} e^{-t} dt \] es la función gamma.
El parámetro \(\alpha\) se denomina parámetro de forma porque es el que determina el grado de asimetría y apuntamiento de la densidad (por ejemplo, el coeficiente de asimetría de una \(\Gamma(\alpha,\beta)\) es \(2/\sqrt{\alpha}\)). El parámetro \(\beta\) se denomina parámetro de escala porque donde más influye es en la dispersión de la distribución. Concretamente, si \(X\sim \Gamma(\alpha,\beta)\) se cumple que \(\mathbb E(X) = \alpha\beta\) y \(\mathbb V(X)=\alpha\beta^2\).
La distribución gamma se utiliza en diversos contextos: fiabilidad, análisis de supervivencia, gestión de riesgos, hidrología, neurología, … Sirve, por ejemplo, para modelizar el coste económico que supone una reclamación para una compañía de seguros, el volumen de precipitación (de lluvia) mensual en una localización, el tiempo de vida útil de un componente o de una maquinaria, …
Comprobar (no necesariamente durante esta sesión práctica) que la función de logverosimilitud de la \(\Gamma(\alpha,\beta)\) es \[ \log L(\alpha,\beta) = -n\log \Gamma(\alpha)-n\alpha\log\beta+(\alpha-1)\sum_{i=1}^n \log x_i -\frac{1}{\beta} \sum_{i=1}^n x_i, \] y que las ecuaciones de verosimilitud son \[\begin{align} \tag{1} 0 & = \frac{\partial}{\partial\alpha} \log L(\alpha,\beta) = -n \frac{\Gamma'(\alpha)}{\Gamma(\alpha)} - n \log\beta + \sum_{i=1}^n \log x_i \label{EcVerAlpha} \\ \tag{2} 0 & = \frac{\partial}{\partial\beta} \log L(\alpha,\beta) = \frac{1}{\beta} \left( -n\alpha+\frac{1}{\beta}\sum_{i=1}^n x_i \right). \label{EcVerBeta} \end{align}\] Es evidente que podemos despejar \(\beta\) de la segunda ecuación obteniendo: \[\begin{equation} \tag{3} \beta = \frac{\bar x}{\alpha}, \label{emvbeta} \end{equation}\] donde \(\bar x=\sum_{i=1}^n x_i /n\) es la media muestral. Si sustituimos la expresión (3) de \(\beta\) en la ecuación de verosimilitud (1) obtenemos la ecuación \[\begin{equation} \tag{4} 0 = - n \psi(\alpha) - n \log \bar x + n\log \alpha + \sum_{i=1}^n \log x_i, \end{equation}\] donde \[\begin{equation} \psi(\alpha) := \frac{d}{d\alpha} \log \Gamma(\alpha) = \frac{\Gamma'(\alpha)}{\Gamma(\alpha)} = \int_0^\infty \left( \frac{e^{-t}}{t} - \frac{e^{-\alpha t}}{1-e^{-t}} \right) dt \end{equation}\] es la llamada función digamma.
Observemos que es imposible obtener una expresión explícita para \(\hat\alpha\), el e.m.v. de \(\alpha\), que es raíz de la ecuación (4). Pero sí podemos aproximar el valor de \(\hat\alpha\) por métodos numéricos. Denotemos \[ g(\alpha) = - n \psi(\alpha) - n \log \bar x + n\log \alpha + \sum_{i=1}^n \log x_i, \] de manera que la ecuación (4) es \(0=g(\alpha)\). El método de Newton aproxima la raíz de esta ecuación mediante la sucesión \[\begin{equation} \label{SucesionNewton1} \tag{5} \alpha_{m} = \alpha_{m-1} - \frac{g(\alpha_{m-1})}{g'(\alpha_{m-1})} = \alpha_{m-1} \left( 1+\frac{\psi(\alpha_{m-1}) + \log\bar x - \log(\alpha_{m-1}) - \overline{\log x}}{1-\alpha_{m-1}\psi'(\alpha_{m-1})} \right), \end{equation}\] donde \(\overline{\log x} := \sum_{i=1}^n \log x_i /n\).
La función \(\psi'\) es la función trigamma. Tanto la función digamma \(\psi\) como la trigamma \(\psi'\) están implementadas en R. Se puede tomar como valor inicial de la iteración el estimador de \(\alpha\) por el método de los momentos \(\alpha_0 = \bar x^2/s^2\).
En la función GammaNewton.R se ha programado la sucesión (5).
GammaNewton = function(muestra,alpha0){
# Sucesion que converge al e.m.v. del parametro alpha de la gamma
m = mean(muestra)
mlog = mean(log(muestra))
alpha1 = alpha0 *(1+(digamma(alpha0)+log(m/alpha0)-mlog)/(1-alpha0*trigamma(alpha0)))
}La función emvalphaGamma.R calcula el e.m.v. del parámetro \(\alpha\).
emvalphaGamma = function(muestra,iter=FALSE){
# Calcula el e.m.v. del parametro alpha de una Gamma a partir de una muestra.
# Si iter=TRUE, imprime en pantalla el error en cada iteracion.
m = mean(muestra)
v = var(muestra)
# Inicio de la iteracion = estimador de alpha por el metodo de momentos:
alpha0 = m^2/v
error = 1
i = 0
while(error>1E-4){
i = i+1
alpha1 = GammaNewton(muestra,alpha0)
error = abs(alpha1-alpha0)
alpha0 = alpha1
if (iter==TRUE){writeLines(paste("Iteracion ",i," Error =",error))}
}
alpha1
}Se puede probar el procedimiento ejecutando el siguiente código:
source("GammaNewton.R")
source("emvalphaGamma.R")
n = 50
alpha = 5
beta = 2
X = rgamma(n,shape=alpha,scale=beta)
m = mean(X)
emvalpha = emvalphaGamma(X,iter=TRUE)
emvbeta = m/emvalpha
writeLines(paste("n =",n))
writeLines(paste("alpha =",alpha," e.m.v.(alpha) =",round(emvalpha,digits=4)))
writeLines(paste("beta =",beta," e.m.v.(beta) =",round(emvbeta,digits=4)))Ejercicio 1: En el fichero WilksTableA2.txt aparecen las precipitaciones medias durante el mes de enero en Ithaca (EEUU) entre los años 1933 y 1981 (Fuente de los datos: Wilks, D.S. (2006). Statistical Methods in the Atmospheric Sciences. Academic Press). La precipitación de lluvia se puede medir en unidades de longitud; concretamente se mide la altura de la lámina de agua recogida en una superficie plana. En este caso las unidades son pulgadas.
Elimina todas las variables del entorno de trabajo con
rm(list=ls())y vuelve a cargar las funciones GammaNewton y emvalphaGamma
source("GammaNewton.R")
source("emvalphaGamma.R")Carga la tabla de datos contenida en WilksTableA2.txt:
Datos = read.table("WilksTableA2.txt",header=FALSE)La primera columna indica el año y la segunda la precipitación media en el correspondiente mes de enero. Extrae la segunda columna
X = Datos$V2y ajústale una distribución gamma por el método de máxima verosimilitud. Dibuja un histograma de los datos y superponle la densidad gamma con valores de los parámetros igual a los emv.
x = seq(0,7,0.01)
hist(X,freq=FALSE,breaks=seq(0,7,0.5))
d = dgamma(x,shape = emvalpha, scale = emvbeta)
lines(x,d,type="l")¿Te parece bueno el ajuste de los datos a la gamma?
La estimación de los e.m.v. de la gamma tiene sesgo (positivo para \(\alpha\) y negativo para \(\beta\)) en la estimación de los parámetros. Vamos a programar un experimento Monte Carlo en el que se extraen 1000 muestras de tamaño \(n\) de una gamma. Para cada muestra Monte Carlo se determinan los e.m.v., \(\hat\alpha\) y \(\hat\beta\), de los parámetros \(\alpha\) y \(\beta\) de la gamma. El resultado final son los promedios de estos estimadores. Comprueba que \(\hat\alpha\) tiende a sobreestimar \(\alpha\) y, por el contrario, \(\hat\beta\) infraestima \(\beta\).
# Experimento Monte Carlo para estudiar el comportamiento
# de los emv de los parametros de la Gamma
n = 50 # Tamaño de cada muestra
alpha = 5
beta = 2
nMC = 1000 # Numero de muestras Monte Carlo
emvalphavector = rep(0,nMC)
emvbetavector = rep(0,nMC)
for (i in 1:nMC){
X = rgamma(n,shape=alpha,scale=beta)
m = mean(X)
emvalpha = emvalphaGamma(X)
emvalphavector[i] = emvalpha
emvbetavector[i] = m/emvalpha
}
writeLines(paste('Tamaño muestral n =',n))
writeLines(paste('Numero de muestras Monte Carlo nMC =',nMC))
writeLines(paste('alpha =',alpha,' e.m.v.(alpha) -> Media =',round(mean(emvalphavector),digits=4),
' Varianza =',round(var(emvalphavector),digits=4)))
writeLines(paste('beta =',beta,' e.m.v.(beta) -> Media =',round(mean(emvbetavector),digits=4),
' Varianza =',round(var(emvbetavector),digits=4)))