Notas para contenedor de docker:
Comando de docker para ejecución de la nota de forma local:
nota: cambiar <ruta a mi directorio> por la ruta de directorio que se desea mapear a /datos dentro del contenedor de docker.
docker run --rm -v <ruta a mi directorio>:/datos --name jupyterlab_r_kernel_local -p 8888:8888 -d palmoreck/jupyterlab_r_kernel:1.1.0password para jupyterlab: qwerty
Detener el contenedor de docker:
docker stop jupyterlab_r_kernel_localDocumentación de la imagen de docker palmoreck/jupyterlab_r_kernel:1.1.0 en liga.
Nota generada a partir de liga
In [1]:
install.packages(c("latex2exp","glmnet"),lib="/usr/local/lib/R/site-library/",
repos="https://cran.itam.mx/")
En esta nota revisamos a la regresión logística con y sin regularización. La regularización que utilizamos es la de Tikhonov también nombrada $\ell_2$ o ridge y la $\ell_1$ o también conocida como lasso (least absolute shrinkage and selection operator, Tibshirani, 1996). Se muestra el uso de métodos de descenso (ver 4.2.Algoritmos_para_UCO) para resolver los problemas de optimización que surgen en los modelos anteriores y no se tiene por objetivo la interpretación de los coeficientes estimados. Se comparan los resultados del paquete glm, glmnet stanford, glmnet cran de R con los obtenidos en la implementación hecha por el prof en algoritmos/R, en específico algoritmos/R/algorithms_for_uco.R para problemas tipo UCO (Unconstrained Convex Optimization).
Sean $\mathcal{C}_0 , \mathcal{C}_1$ dos clases ajenas y $x \in \mathbb{R}^n$. El problema de clasificación consiste en clasificar al vector $x$ en alguna de las dos clases anteriores de modo que se minimice el error de clasificación.
Ejemplos de lo anterior los encontramos en medicina (persona enferma o no dada una serie de mediciones en sangre), finanzas (persona sujeta a un crédito bancario o no dado un historial crediticio) o clasificación de textos (spam o no spam).
En esta nota revisamos a la regresión logística que define un modelo de probabilidad de pertenenecia a una clase.
In [1]:
library(ggplot2)
library(latex2exp)
El modelo por regresión logística tiene por objetivo modelar las probabilidades de pertenencia a cada una de las clases $\mathcal{C}_0, \mathcal{C}_1$ dado el vector de atributos $ x \in \mathbb{R}^n$: $p(\mathcal{C}_0|x) , p(\mathcal{C}_1|x)$.
En la regresión logística se utiliza la función sigmoide $\sigma:\mathbb{R} \rightarrow \mathbb{R}$:
$$\sigma(t)=\frac{1}{1+\exp(-t)}$$para modelar ambas probabilidades ya que mapea todo el eje real al intervalo $[0,1]$. Además resulta ser una aproximación continua y diferenciable a la función de Heaviside $H:\mathbb{R} \rightarrow \mathbb{R}$
$$H(t) = \begin{cases} 1 & \text{si } t \geq 0,\\ 0 & \text{si } t <0\\ \end{cases} $$
In [2]:
t <- seq(from = -10, to = 10, by = .01)
In [3]:
Heaviside <- function(t) 1*(t>0)
In [4]:
options(repr.plot.width=4, repr.plot.height=4) #esta línea sólo se ejecuta para jupyterlab con R
In [5]:
qplot(t, Heaviside(t))
In [6]:
options(repr.plot.width=7.5, repr.plot.height=7.5) #esta línea sólo se ejecuta para jupyterlab con R
A continuación graficamos a la sigmoide $\sigma(ht)$ para distintos valores de $h \in \{-3, -1, -1/2, 1/2, 1, 3\}$:
In [7]:
sigmoide <- function(t) 1/(1+exp(-t))
In [8]:
h<-c(-3,-1,-1/2,1/2,1,3)
h_reps <- as.vector(
vapply(1:length(h),function(i) rep(h[i], NROW(t)),
numeric(length(t))
)
)
sigmoide_reps <- as.vector(
vapply(1:length(h), function(i) sigmoide(h[i]*t),
numeric(length(t))
)
)
In [9]:
df <- data.frame(x = rep(t, NROW(h)),
y = sigmoide_reps,
h = as.factor(h_reps)
)
In [10]:
gg <- ggplot(df)
In [11]:
gg + geom_line(aes(x = x, y = y, group = h, color = h)) +
xlab('t') +
ggtitle(TeX('$\\sigma(ht)$'))
Obsérvese la forma de cada curva al variar $h$ en la función $\sigma(ht)$. Una regla de clasificación podría ser clasificar como perteneciente a $\mathcal{C}_0$ si la probabilidad (modelada por la curva sigmoide) es menor a $0.25$ (punto de corte) y perteneciente a $\mathcal{C}_1$ si es mayor o igual a $0.25$. Para diferentes curvas sigmoides presentadas en la gráfica anterior obsérvese que al fijar el punto de corte y tomar un valor de $t$ en el eje horizontal, la pertenencia a alguna de las clases es menos sensible al variar $t$ que en otras curvas.
Así, la función sigmoide permite modelar la probabilidad de pertenencia a la clase $\mathcal{C}_1:$
$$p(\mathcal{C}_1| x)=\sigma(a)$$para alguna $a \in \mathbb{R}$.
Con el teorema de Bayes se obtiene el valor de $a$:
$$ \begin{eqnarray} p(\mathcal{C}_1|x) &=& \frac{p(x|\mathcal{C}_1)p(\mathcal{C}_1)}{p(x|\mathcal{C}_0)p(\mathcal{C}_0)+p(x|\mathcal{C}_1)p(\mathcal{C}_1)} \nonumber \\ &=& \left ( 1+ \frac{p(x|\mathcal{C}_0)p(\mathcal{C}_0)}{p(x|\mathcal{C}_1)p(\mathcal{C}_1)} \right )^{-1} \nonumber \end{eqnarray} $$Por lo tanto:
Comentarios:
se le conoce como logit y modela el log momio:
que tiene una interpretación directa en términos de las probabilidades de pertenencia a cada clase $\mathcal{C}_0,\mathcal{C}_1$.
De forma similar como en el modelo por mínimos cuadrados lineales (ver 4.3.Minimos_cuadrados_R) se modeló a la variable respuesta $y$ con una función lineal en sus parámetros, en el modelo en regresión logística con dos clases se propone una función lineal en un vector de parámetros $(\beta_0,\beta) \in \mathbb{R}^{n+1}$ definida por el logit:
Obsérvese que si $y$ es considerada como variable respuesta que está en función de $x \in \mathbb{R}^{n+1}$ dado el vector $(\beta_0, \beta)$ se tiene:
que se lee "la probabilidad de pertenencia a la clase $\mathcal{C}_1$ dado el vector de atributos $x$ es igual a $y$".
Comentarios:
con $x \in \mathbb{R}$.
La notación $y(x | \beta_0, \beta)$ se utiliza para denotar que $(\beta_0, \beta)$ es un vector de parámetros a estimar, en específico $\beta_0, \beta_1, \dots, \beta_n$, esto es: $n+1$ parámetros a estimar.
La variable de optimización es $(\beta_0, \beta) \in \mathbb{R}^{n+1}$.
Dados $(x_0,\hat{y}_0), \dots (x_m, \hat{y}_m)$ puntos se desean modelar $m+1$ probabilidades de pertenencias a las clases $\mathcal{C}_0, \mathcal{C}_1$ representadas con las etiquetas $\hat{y}_i \in \{0,1\} \forall i=0,1,\dots, m$. El número $0$ representa a la clase $\mathcal{C}_0$ y el $1$ a la clase $\mathcal{C}_1$. El vector $x_i \in \mathbb{R}^n$ .
Cada probabilidad se modela como $y_0=y_0(x_0|\beta_0, \beta),y_1=y_1(x_1|\beta_0, \beta),\dots,y_m=y_m(x_n|\beta_0, \beta)$ utilizando:
$$p(\mathcal{C}_1|x_i) = y_i(x_i|\beta_0,\beta) = \frac{1}{1+ e^{-(\beta_0 + \beta^T x)}} \forall i=0,1,\dots,m.$$Los $n+1$ parámetros $\beta_0, \beta_1, \dots, \beta_n$ se ajustan maximizando la función de verosimilitud:
$$ \mathcal{L}(\beta_0, \beta|x)=\displaystyle \prod_{i=0}^n y_i^{\hat{y}_i}(1-y_i)^{1-\hat{y}_i} $$donde: $\hat{y}_i \sim \text{Bernoulli}(y_i)$ y por tanto $\hat{y}_i \in \{0,1\}$: $\hat{y}_i = 1$ con probabilidad $y_i$ y $\hat{y}_i = 0$ con probabilidad $1-y_i$. Entonces se tiene el problema:
Lo anterior es equivalente a maximizar la log-verosimilitud:
$$ \begin{eqnarray} \ell(\beta_0, \beta |x)&=&\log(\mathcal{L}(\beta_0, \beta| x)) \nonumber\\ &=&\displaystyle \sum_{i=1}^m \hat{y}_i\log(y_i)+(1-\hat{y}_i)\log(1-y_i) \nonumber\\ &=&\displaystyle \sum_{i=1}^m\hat{y}_i (\beta_0, \beta)^T x_i-\log(1+\exp((\beta_0, \beta)^Tx_i) \nonumber \end{eqnarray} $$o a minimizar la devianza:
Comentario:
La devianza es una función convexa pues su Hessiana* es:
$$ \nabla^2 D(\beta_0, \beta |x) = 2A^TPA $$con: $P$ una matriz diagonal con entradas $y_i(1-y_i)$ donde: $y_i$ está en función de $x_i$: $ y_i(x_i|\beta_0,\beta) = \frac{1}{1+ e^{-(\beta_0 + \beta^T x_i)}} \forall i=0,1,\dots,m $ y la matriz A es:
*La expresión anterior de la Hessiana se obtiene a partir de la expresión del gradiente:
$$ \nabla D(\beta_0, \beta|x) = 2 \displaystyle \sum_{i=1}^m \left( y_i - \hat{y}_i \right )x_i = 2\sum_{i=1}^m \left( p(\mathcal{C}_1|x_i) - \hat{y}_i \right )x_i = 2A^T(p-\hat{y}) $$donde:
Así, la Hessiana de la devianza es simétrica semidefinida positiva y por tanto es una función convexa. Ver el apéndice de definiciones para funciones convexas de 4.1.Optimizacion_numerica_y_machine_learning.
Por lo anterior podemos utilizar los métodos de descenso para resolver problemas de optimización convexa sin restricciones revisados en la nota 4.2.Algoritmos_para_UCO.
Como primer ejemplo utilizamos el conocido dataset de Iris en el que se muestran tres especies del género Iris. Las especies son: I. setosa, I. virginica y I. versicolor:
Imagen obtenida de Iris Dataset.
In [12]:
library(datasets)
library(dplyr)
library(glmnet)
In [13]:
data(iris)
summary(iris)
In [14]:
print("Número de renglones y columnas:")
nrow(iris)
ncol(iris)
Para el modelo de regresión logística en dos clases utilizamos a la clase virginica y versicolor.
In [15]:
iris_subset <- filter(iris, Species %in% c("virginica", "versicolor"))
In [16]:
head(iris_subset)
In [17]:
nrow(iris_subset)
Obsérvese que la variable Sepal.Length y Petal.Length están altamente correlacionadas por lo que quitaremos a Petal.length:
In [21]:
cor(as.matrix(iris_subset[,-ncol(iris_subset)]))
In [22]:
A<-iris_subset[,-c(3,ncol(iris_subset))]
In [23]:
head(A)
La variable $\hat{y}$ estará dada por las clases versicolor y virginica:
In [24]:
head(iris_subset[,ncol(iris_subset)])
In [25]:
tail(iris_subset[,ncol(iris_subset)])
Transformamos a tipo numérica la columna Species:
In [26]:
iris_subset$Clase <- as.numeric(iris_subset$Species %in% 'versicolor')
In [27]:
head(iris_subset)
In [28]:
tail(iris_subset)
La clase $\mathcal{C}_1$ es versicolor y $\mathcal{C}_0$ es virginica.
In [29]:
glm.out <- glm(Clase ~ Sepal.Length + Sepal.Width + Petal.Width,
data = iris_subset,family = binomial)
In [30]:
summary(glm.out)
Y nos interesan las estimaciones de los coeficientes $\beta_0, \beta_1, \beta_2, \beta_3$:
In [31]:
print(glm.out$coefficients)
La función objetivo como se revisó en la sección anterior está dada por la expresión de la devianza:
donde: $\hat{y}_i \in \{0,1\}$, $x_i$ $i$-ésimo renglón de matriz $A \in \mathbb{R}^{100 \times 3}$.
Añadimos la columna que indica uso de intercepto y por tanto de un modelo con $\beta_0$:
In [32]:
A<-cbind(rep(1,nrow(A)),A)
In [33]:
head(A)
In [34]:
y_hat <- iris_subset$Clase
La clase $\mathcal{C}_1$ es versicolor:
In [35]:
head(y_hat)
La clase $\mathcal{C}_0$ es virginica:
In [36]:
tail(y_hat)
Cargamos las funciones programadas en el directorio algoritmos/R:
In [37]:
#load numerical differentiation
#load utils
#load algorithms for unconstrained convex optimization
#load line search
dir_R="algoritmos/R"
source(paste(dir_R,"/numerical_differentiation.R",
sep=""))
source(paste(dir_R,"/utils.R",
sep=""))
source(paste(dir_R,"/algorithms_for_uco.R",
sep=""))
source(paste(dir_R,"/line_search.R",
sep=""))
Convertimos a tipo matriz:
In [38]:
A<-as.matrix(A)
In [39]:
head(A)
Función objetivo:
In [40]:
fo <- function(beta){
mat_vec <- A%*%beta
2*sum(log(1+exp(mat_vec))-y_hat*mat_vec)
}
Punto inicial $\beta ^{(0)}:$
In [41]:
beta_0 <- c(1,1,1,1)
$\beta^*$:
In [42]:
beta_ast_glm <- glm.out$coefficients
In [43]:
print(beta_ast_glm)
In [44]:
p_ast <- fo(beta_ast_glm)
In [45]:
p_ast
Argumentos para el método de Newton:
In [46]:
tol <- 1e-8
tol_backtracking <- 1e-14
maxiter <- 30
In [47]:
l<-Newtons_method(fo, beta_0, tol, tol_backtracking, beta_ast_glm, p_ast, maxiter)
In [48]:
beta <- l[[1]]
total_of_iterations <- l[[2]]
Err_plot <- l[[3]]
beta_plot <- l[[4]]
Secuencia de minimización $\beta^{(k)}$:
In [49]:
beta_plot
$\beta$ aproximada por el método de Newton:
In [50]:
print(beta)
$\beta^*$:
In [51]:
print(beta_ast_glm)
In [52]:
compute_error(x_obj = beta_ast_glm,x_approx = beta)
Tenemos alrededor de 3 dígitos de precisión.
Comparemos con la función glmnet del paquete del mismo nombre (no usamos intercepto):
In [54]:
fit<-glmnet(A[,-1],y_hat,family="binomial",alpha=0,lambda=0,standardize=F,nlambda=1, thresh=1e-8)
In [55]:
beta_ast_glmnet <- as.matrix(fit$beta)
In [56]:
print(beta_ast_glmnet)
Error relativo:
In [57]:
compute_error(x_obj = beta_ast_glmnet,x_approx = beta[-1])
Tenemos alrededor de 3 dígitos de precisión.
Para individuo $i$ se tiene:
Por ejemplo para el renglón $1$ de $A$: (que sabemos que pertenece a la clase $\mathcal{C}_1$: versicolor):
Estimación de la probabilidad de pertenencia a la clase $\mathcal{C}_1$: versicolor con lo estimado por glm:
In [58]:
1/(1+exp(-sum(A[1,]*beta_ast_glm)))
Estimación de la probabilidad de pertenencia a la clase $\mathcal{C}_0$: virginica con lo estimado por glm:
In [59]:
exp(-sum(A[1,]*beta_ast_glm))/(1+exp(-sum(A[1,]*beta_ast_glm)))
Estimación de la probabilidad de pertenencia a la clase $\mathcal{C}_1$:versicolor con lo estimado por descenso por Newton:
In [60]:
1/(1+exp(-sum(A[1,]*beta)))
Estimación de la probabilidad de pertenencia a la clase $\mathcal{C}_0$:virginica con lo estimado por descenso por Newton:
In [61]:
exp(-sum(A[1,]*beta))/(1+exp(-sum(A[1,]*beta)))
Por ejemplo para el último renglón de $A$: (que sabemos que pertenece a la clase $\mathcal{C}_0$: virginica):
Estimación de la probabilidad de pertenencia a la clase $\mathcal{C}_1$: versicolor con lo estimado por glmn:
In [62]:
1/(1+exp(-sum(A[nrow(A),]*beta_ast_glm)))
Estimación de la probabilidad de pertenencia a la clase $\mathcal{C}_0$: virginica con lo estimado por glmn:
In [63]:
exp(-sum(A[nrow(A),]*beta_ast_glm))/(1+exp(-sum(A[nrow(A),]*beta_ast_glm)))
Estimación de la probabilidad de pertenencia a la clase $\mathcal{C}_1$:versicolor con lo estimado por descenso por Newton:
In [64]:
1/(1+exp(-sum(A[nrow(A),]*beta)))
Estimación de la probabilidad de pertenencia a la clase $\mathcal{C}_0$:virginica con lo estimado por descenso por Newton:
In [65]:
exp(-sum(A[nrow(A),]*beta))/(1+exp(-sum(A[nrow(A),]*beta)))
A continuación consideramos regularización con ridge y lasso, ver 4.3.Minimos_cuadrados_R.ipynb.
In [66]:
reg<-0.5
In [67]:
fit<-glmnet(A,y_hat,family="binomial",alpha=0,lambda=reg,standardize=F,nlambda=1, thresh=1e-8)
$\beta^*:$
In [68]:
beta_ast<-as.matrix(fit$beta)
In [69]:
print(beta_ast)
Función objetivo de la función glmnet con regularización ridge:
In [70]:
mpoints<-nrow(A)
In [72]:
fo <-function(beta){
mat_vec <- A %*% beta
1/mpoints*(sum(log(1 + exp(mat_vec)) - y_hat * mat_vec)) + reg/2*sum(beta[2:length(beta)]*beta[2:length(beta)])
}
Punto inicial $\beta ^{(0)}:$
In [73]:
beta_0 <- c(0,0,0,0)
In [74]:
p_ast <- fo(beta_ast)
In [75]:
p_ast
In [76]:
tol <- 1e-12
tol_backtracking <- 1e-14
maxiter <- 30
In [77]:
l<-Newtons_method(fo, beta_0, tol, tol_backtracking, beta_ast, p_ast, maxiter)
In [78]:
beta <- l[[1]]
total_of_iterations <- l[[2]]
Err_plot <- l[[3]]
beta_plot <- l[[4]]
$\beta$ aproximada por el método de Newton:
In [79]:
print(beta)
In [80]:
print(beta_ast)
Error relativo:
In [81]:
compute_error(beta_ast[-1],beta[-1])
Sin considerar el intercepto tenemos alrededor de $6$ dígitos de precisión.
In [82]:
A<-A[,-1]
In [83]:
head(A)
Solución con glmnet:
In [84]:
fit<-glmnet(A,y_hat,family="binomial",alpha=0,lambda=reg,standardize=F,nlambda=1, intercept=F, thresh=1e-8)
In [85]:
beta_ast <- as.matrix(fit$beta)
$\beta^*:$
In [86]:
print(beta_ast)
Solución vía método de Newton:
Punto inicial $\beta ^{(0)}:$
In [87]:
beta_0 <- c(0,0,0)
In [88]:
p_ast <- fo(beta_ast)
In [89]:
p_ast
In [90]:
tol <- 1e-12
tol_backtracking <- 1e-14
maxiter <- 30
In [91]:
fo <-function(beta){
mat_vec <- A %*% beta
1/mpoints*(sum(log(1 + exp(mat_vec)) - y_hat * mat_vec)) + reg/2*sum(beta*beta)
}
In [92]:
l<-Newtons_method(fo, beta_0, tol, tol_backtracking, beta_ast, p_ast, maxiter)
In [93]:
beta <- l[[1]]
total_of_iterations <- l[[2]]
Err_plot <- l[[3]]
beta_plot <- l[[4]]
$\beta$ aproximada por el método de Newton:
In [94]:
print(beta)
In [95]:
print(beta_ast)
Error relativo:
In [96]:
compute_error(x_approx = beta, x_obj = beta_ast)
Tenemos alrededor de $5$ dígitos de precisión.
Función objetivo de la función glmnet con regularización lasso:
In [97]:
reg <- 0.01
Solución vía glmnet:
In [98]:
fit<-glmnet(A,y_hat,family="binomial",alpha=1,lambda=reg,standardize=F,nlambda=1, intercept=F, thresh=1e-8)
In [99]:
beta_ast <- as.matrix(fit$beta)
$\beta^*:$
In [100]:
print(beta_ast)
Solución vía método de Newton:
Usaremos la siguiente función quita_signo que ayuda a aproximar la derivada de la función objetivo $f_o$ y así lidiar con la no diferenciabilidad en $0$.
La función quita_signo realiza:
obsérvese que el valor elegido en el tercer caso es el más pequeño positivo normalizado en un sistema de punto flotante, ver 1.2.Sistema_de_punto_flotante.
In [101]:
quita_signo<-function(beta){
beta<-sign(beta)*beta
ind <- beta < .Machine$double.xmin & beta > -.Machine$double.xmin
beta[ind] <- .Machine$double.xmin
beta
}
Comentario: la definición de la función quita_signo se basa en lo que se conoce como subdifferential que es un conjunto de subderivatives, útiles para generalizar las derivadas para funciones convexas que no son diferenciables en puntos de su dominio.
In [103]:
fo <-function(beta){
mat_vec <- A %*% beta
1/mpoints*(sum(log(1 + exp(mat_vec)) - y_hat * mat_vec)) + reg*sum(quita_signo(beta))
}
Punto inicial $\beta ^{(0)}:$
In [104]:
beta_0<-c(1,1,1)
In [105]:
l<-Newtons_method(fo, beta_0, tol, tol_backtracking, beta_ast, p_ast, maxiter)
In [106]:
beta <- l[[1]]
total_of_iterations <- l[[2]]
Err_plot <- l[[3]]
beta_plot <- l[[4]]
$\beta$ aproximada por el método de Newton:
In [107]:
print(beta)
In [108]:
print(beta_ast)
Error relativo:
In [109]:
compute_error(x_approx = beta, x_obj = beta_ast)
Tenemos alrededor de 3 dígitos de precisión.