Search

I’m your maximum likelihood estimate

En mi anterior post,  Descubriendo la regresión logística, dejaba indicado que, para estimar los coeficientes del modelo de regresión logística, se aplica el método de máxima verosimilitud (maximum likelihood, en inglés). En realidad este método se puede utilizar para calcular los parámetros de cualquier otro modelo.

¿En qué consiste este método?

Bueno, para una muestra aleatoria de sucesos independientes, idénticamente distribuidos (x1,…,xn), cuya distribución depende de un parámetro θ, la función de verosimilitud se define así:

    \[ \mathcal{L}(\theta ) = f(x_{1};\theta)\cdots f(x_{n};\theta) =  \prod_{i=1}^{n}f(x_{i};\theta) \]

Por conveniencia, para hacer más sencillos los cálculos, se suele utilizar el logaritmo de la función de verosimilitud, también llamada (log)verosimilitud, que tiene esta pinta:

    \[ \ell(\theta) = \ln\mathcal{L}(\theta) = \ln\prod_{i=1}^{n}f(x_{i};\theta) = \sum_{i=1}^{n}\ln f(x_{i};\theta) \]

Con el método de máxima verosimilitud, obtendremos el valor del parámetro θ que maximice la función de verosimilitud anterior. Al parámetro θ se le llama estimación de máxima verosimilitud, (EMV).

Como ya sabemos, en la regresión logística la variable respuesta, yi, sigue una distribución de Bernoulli tal que:

    \[ y_{i} = \begin{cases} 0 & P_{i}(0) = 1-p_{i}\\ 1 & P_{i}(1) = p_{i} \end{cases} \]

y su función de probabilidad tiene esta pinta:

    \[ f(y_{i}) = p_{i}^{y_{i}} \cdot (1-p)_{i}^{1-y_{i}} \]

Si recordamos, dos de las expresiones que obtuvimos en Descubriendo la regresión logística eran:

(1)   \begin{equation*}  \ln\begin{pmatrix}\frac{p_{i}}{1-p_{i}}\end{pmatrix}= B\cdot x_{i} =b_{0}\cdot x_{i0} + b_{1}\cdot x_{i1}+\cdots + b_{p}\cdot x_{ip} \end{equation*}

(2)   \begin{equation*}  p_{i} = \frac{1}{1 + e^{-(B\cdot x_{i})}} \end{equation*}

donde pi es la probabilidad de que ocurra el evento, xip son los factores y bi los parámetros asociados a cada factor, que, en principio, están relacionados con la probabilidad de que el suceso ocurra. Por ejemplo, para dos únicos parámetros tales que b0 = a y xi0 = 1  y b1 = b  tendríamos la expresión para un regresión lineal simple:

    \[\ln\begin{pmatrix}\frac{p_{i}}{1-p_{i}}\end{pmatrix}=a + b\cdot x_{i}\]

Así pues, B son los parámetros que necesitamos estimar.

Entonces, partiendo de la definición para la función de verosimilitud hacemos lo siguiente:

    \[ \mathcal{L}(b) = \prod_{i=1}^{n}f(y_{i})=\prod_{i=1}^{n}p_{i}^{y_{i}} \cdot (1-p)_{i}^{1-y_{i}}\Rightarrow \]

    \[ \rightarrow \enspace Aplicamos \enspace logaritmos \]

    \[ \ln\mathcal{L}(b)  = \ell(b) = \ln \prod_{i=1}^{n}p_{i}^{y_{i}} \cdot (1-p)_{i}^{1-y_{i}}= \sum_{i=1}^{n} y_{i}\ln p_{i} + (1-y_{i})\ln(1-p_{i})\Rightarrow \]

    \[\Rightarrow \ell(b) =\sum_{i=1}^{n} y_{i}\ln p_{i}+\ln(1-p_{i})-y_{i}\ln(1-p_{i}) = \sum_{i=1}^{n}y_{i}(\ln p_{i}-\ln(1-p_{i}) ) + \ln(1-p_{i})\Rightarrow\]

    \[ \begin{array}{lcl} \rightarrow \enspace Aplicamos \enspace que \enspace  \ln p_{i}-\ln(1-p_{i}) = \ln\begin{pmatrix}\frac{p_{i}}{1-p_{i}}\end{pmatrix}  \\ \rightarrow  \enspace Aplicamos \enspace la \enspace ecuaci\acute{o}n \enspace  (1) \]

    \[\Rightarrow \ell(b) = \sum_{i=1}^{n}y_{i}\ln\begin{pmatrix}\frac{p_{i}}{1-p_{i}}\end{pmatrix} + \ln(1-p_{i}) = \sum_{i=1}^{n}y_{i}(B\cdot x_{i}) + \ln(1-p_{i}) \Rightarrow\]

    \[ \begin{array}{lcl} \rightarrow \enspace Aplicamos \enspace lo \enspace  siguiente, utilizando \enspace la \enspace ecuaci\acute{o}n \enspace  (2): \\ 1-p_{i}=1-\frac{1}{1 + e^{-(B\cdot x_{i})}}=\frac{1-e^{-(B\cdot x_{i})}+1}{1 + e^{-(B\cdot x_{i}}}=\frac{e^{-(B\cdot x_{i})}}{1+e^{-(B\cdot x_{i})}}= \frac{1}{e^{(B\cdot x_{i})}+1} \]

    \[\Rightarrow \ell(b) = \sum_{i=1}^{n}y_{i}\cdot B\cdot x_{i}+ \ln\begin{pmatrix}\frac{1}{e^{(B\cdot x_{i})}+1}\end{pmatrix})=\sum_{i=1}^{n}y_{i}\cdot B\cdot x_{i}+\ln1 - \ln(e^{(B\cdot x_{i})}+1) \Rightarrow\]

    \[ \Rightarrow \boxed{\ell(b) = \sum_{i=1}^{n}y_{i}\cdot B\cdot x_{i}- \ln(e^{(B\cdot x_{i})}+1)} \]

Una vez que tenemos nuestra función (log)verosimilitud, tenemos que maximizarla.

Para ello lo primero que tenemos que hacer es calcular las derivadas parciales de la función respecto a los parámetros B:

    \[\frac{\partial \ell(b)}{\partial b}=\frac{\partial \sum_{i=1}^{n}y_{i}\cdot B\cdot x_{i} - \ln(e^{(B\cdot x_{i})}+1)}{\partial b}=\sum_{i=1}^{n} y_{i}\cdot x_{i} - \frac{e^{(B\cdot x_{i})}}{e^{(B\cdot x_{i})}+1}= \sum_{i=1}^{n} y_{i}\cdot x_{i} - \frac{1}{e^{-(B\cdot x_{i})}+1} \Rightarrow\]

    \[ \begin{array}{lcl} \rightarrow Aplicamos: \\ p_{i} = \frac{1}{1 + e^{-(B\cdot x_{i})}} \]

    \[\Rightarrow \frac{\partial \ell(b)}{\partial b} = \sum_{i=1}^{n} y_{i}\cdot x_{i} - p_{i}\cdot x_{i} = \]

La expresión para la primera derivada podemos escribirla en notación matricial así:

    \[ \boxed{\frac{\partial \ell(b)}{\partial b} = X(y - p)} \]

Ahora vamos a calcular la segunda derivada de la función:

    \[\frac{\partial^2 \ell(b)}{\partial b^2} = -x_{i}\cdot \frac{\partial p_{i}}{\partial b} \Rightarrow \]

    \[ \begin{array}{lcl} \rightarrow Calculamos  \enspace \frac{\partial p_{i}}{\partial b}\\ \frac{\partial p_{i}}{\partial b}=\frac{\partial }{\partial b}\left ( \frac{1}{1 + e^{-(B\cdot x_{i})}} \right )= \frac{-e^{-(B\cdot x_{i})}\cdot x_{i}}{\left ( 1 + e^{-(B\cdot x_{i})} \right )^{2}} = -\frac{1}{1 + e^{-(B\cdot x_{i})}}\cdot\frac{e^{-(B\cdot x_{i})}}{1 + e^{-(B\cdot x_{i})}}\cdot x_{i} = -p_{i}\cdot (1-p_{i})\cdot x_{i} \]

    \[\Rightarrow \frac{\partial^2 \ell(b)}{\partial b^2} = \sum_{i=1}^{n}x_{i}\cdot p_{i}\cdot (1-p_{i})\cdot x_{i}\]

En notación matricial sería:

    \[ \boxed{\frac{\partial^2 \ell(b)}{\partial b^2} = X\cdot W \cdot X^{T}} \]

donde W es la matriz diagonal:

    \[ W = \begin{pmatrix} \cdots & & \\ & p_{i}(1-p_{i}) & \\ & & \cdots \end{pmatrix} \]

Para calcular los extremos (máximos y mínimos) de una función, hay que igualar a cero la primera derivada de dicha función, y despejar los coeficientes. Pero, dependiendo de la función, ésto no siempre es sencillo. Para ayudarnos en esa tarea, existe otro método, el método de Newton-Raphson. Con éste método podemos aproximar los valores de las raíces de cualquier función. Es un método iterativo que utiliza esta ecuación para actualizar los parámetros en cada iteración:

    \[b_{n+1}= b_{n} - \frac{f(b_{n})}{{f}'(b_{n})} \Rightarrow\]

    \[ \Rightarrow\Delta b=\frac{f(b_{n})}{{f}'(b_{n})}\]

En nuestro caso, teniendo en cuenta que f(bn) sería igual a la expresión de nuestra primera derivada y f'(bn) sería igual a la expresión de la segunda derivada, la iteración de Newton-Raphson se escribiría así:

    \[ \boxed{\Delta b=(XWX^{T})^{-1}X(y-p)} \]

El algoritmo de estimación de los coeficientes B es el siguiente:

  1. Damos un valor al coeficiente b = 0.
  2. Se estima el vector p de probabilidades para cada evento (que en esta primera iteración serán pi =0.5) y la matriz W.
  3. Se calcula el vector z.
  4. Se actualiza b.

Los pasos 2 al 4 se repiten hasta lograr la convergencia. Aunque ésta no se garantiza siempre, lo normal es alcanzarla.

Bueno, pues vamos a dejar la teoría y vamos a intentar poner todo lo que hemos visto hasta ahora en práctica, desarrollando un programita en Python para hacer todos estos cálculos, y obtener los coeficientes de una regresión logística.

Empezamos importando los módulos que vamos a necesitar:

Lo siguinte que haremos será generar unos datos, que si no ,no hacemos nada :-). Lo haremos con la función genDataset que cree en el post Atrévete a crear tu propio dataset. Vamos a crear un dataset con 4 características y 50000 registros, con valores entre 0 y 10, por ejemplo.

Echamos un vistazo a los datos:

Bueno, ahora que ya tenemos unos datos preparados, vamos a crear el modelo.

Definiremos una clase, a la que llamaremos LogReg_mv:

Añadiremos los siguientes métodos:

El método w_matrix para calcular la matriz W:

El método prob_matrix para calcular la probabilidad asociada a cada suceso:

El método fit para el cálculo de los coeficientes:

Y por último, el método predict, para predecir los valores de la variable respuesta:

Con lo que nuestra clase tendría esta pinta:

El siguiente paso será utilizar la función train_test_split para dividir el dataset en dos datasets, uno de entrenamiento y otro de prueba:

Genial. Ya podemos entrenar a nuestro modelo:

Nos dice que ha necesitado 4 iteraciones para conseguir esos coeficientes.

Podemos sacar una gráfica con el incremento de cada iteración:

Le pedimos al modelo la variable respuesta para el dataset X_test

Evaluamos el modelo. Podemos utilizar la matriz de confusión:

Bueno, con un 73% de accuracy, el modelo es más bien ‘normalito’.

Vamos a utilizar el módulo statmodels para, con los mismos datos, crear un modelo logístico y obtener los coeficientes, para ver si coinciden con los que ha calculado nuestro modelo.

Aunque nuestro modelo puede mejorar clasificando, los coeficientes de la regresión logística si que los calcula bastante bien.

Comparison of the coefficients 
Coef.LogReg_mvstatsmodels
f1-0.179-0.1794
f20.0690.0690
f3-0.1691 -0.1691
f40.06660.0666

¡Anímate y prueba con tu propio dataset!

Puedes descargarte el código de mi 

Deja un comentario

Este sitio usa Akismet para reducir el spam. Aprende cómo se procesan los datos de tus comentarios.