<jupyter_start><jupyter_text>**Notas para contenedor de docker:** Comando de docker para ejecución de la nota de forma local:nota: cambiar `` por la ruta de directorio que se desea mapear a `/datos` dentro del contenedor de docker.```docker run --rm -v :/datos --name jupyterlab_r_kernel_local -p 8888:8888 -d palmoreck/jupyterlab_r_kernel:1.1.0```password para jupyterlab: `<PASSWORD>`Detener el contenedor de docker:```docker stop jupyterlab_r_kernel_local``` Documentación de la imagen de docker `palmoreck/jupyterlab_r_kernel:1.1.0` en [liga](https://github.com/palmoreck/dockerfiles/tree/master/jupyterlab/r_kernel). --- Nota generada a partir de [liga](https://www.dropbox.com/s/6isby5h1e5f2yzs/4.2.Problemas_de_optimizacion_convexa.pdf?dl=0)<jupyter_code>install.packages(c("latex2exp","glmnet"),lib="/usr/local/lib/R/site-library/",
                repos="https://cran.itam.mx/")<jupyter_output>also installing the dependencies ‘iterators’, ‘foreach’, ‘shape’<jupyter_text>Problema de clasificación En esta nota revisamos a la **regresión logística con y sin regularización**. La **regularización** que utilizamos es la de **[Tikhonov](https://en.wikipedia.org/wiki/Tikhonov_regularization)** también nombrada $\ell_2$ o ***ridge*** y la $\ell_1$ o también conocida como **[*lasso*](https://en.wikipedia.org/wiki/Lasso_(statistics)** (*least absolute shrinkage and selection operator*, Tibshirani, 1996). Se muestra el uso de **métodos de descenso** (ver [4.2.Algoritmos_para_UCO](https://github.com/ITAM-DS/analisis-numerico-computo-cientifico/blob/master/temas/IV.optimizacion_convexa_y_machine_learning/4.2.Algoritmos_para_UCO.ipynb)) 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](https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/glm), [glmnet stanford](https://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html), [glmnet cran](https://cran.r-project.org/web/packages/glmnet/index.html) de R con los obtenidos en la implementación hecha por el prof en [algoritmos/R](algoritmos/R), en específico [algoritmos/R/algorithms_for_uco.R](https://github.com/ITAM-DS/analisis-numerico-computo-cientifico/blob/master/temas/IV.optimizacion_convexa_y_machine_learning/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.  Regresión logística: clasificación en $\mathcal{C}_0, \mathcal{C}_1$<jupyter_code>library(ggplot2)
library(latex2exp)<jupyter_output><empty_output><jupyter_text>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](https://en.wikipedia.org/wiki/Sigmoid_function)** $\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](https://en.wikipedia.org/wiki/Heaviside_step_function)** $H:\mathbb{R} \rightarrow \mathbb{R}$$$H(t) = \begin{cases}1 & \text{si } t \geq 0,\\0 & \text{si } t <0\\\end{cases}$$<jupyter_code>t <- seq(from = -10, to = 10, by = .01)
Heaviside <- function(t) 1*(t>0)
options(repr.plot.width=4, repr.plot.height=4) #esta línea sólo se ejecuta para jupyterlab con R
qplot(t, Heaviside(t))
options(repr.plot.width=7.5, repr.plot.height=7.5) #esta línea sólo se ejecuta para jupyterlab con R<jupyter_output><empty_output><jupyter_text>A continuación graficamos a la sigmoide $\sigma(ht)$ para distintos valores de $h \in \{-3, -1, -1/2, 1/2, 1, 3\}$:<jupyter_code>sigmoide <- function(t) 1/(1+exp(-t))
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))      
                                  )
                           )
df <- data.frame(x = rep(t, NROW(h)), 
                 y = sigmoide_reps, 
                 h = as.factor(h_reps)
                 )
gg <- ggplot(df)
gg + geom_line(aes(x = x, y = y, group = h, color = h)) + 
   xlab('t') + 
   ggtitle(TeX('$\\sigma(ht)$'))<jupyter_output><empty_output><jupyter_text>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](https://en.wikipedia.org/wiki/Bayes%27_theorem) 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: $$\begin{eqnarray}a(x)&=&\log\left( \frac{p(x|\mathcal{C}_1)p(\mathcal{C}_1)}{p(x|\mathcal{C}_0)p(\mathcal{C}_0)} \right ) \nonumber\end{eqnarray}$$ **Comentarios:*** Entre algunas propiedades que tiene la función $\sigma(\cdot)$ se encuentran:$$\begin{eqnarray}\sigma (-t)&=&1-\sigma (t) \nonumber \\\frac{d\sigma (t)}{dt}&=&\sigma (t)(1-\sigma (t)) \nonumber\end{eqnarray}$$ * En Estadística a la función:$$a=\log\left(\frac{\sigma}{1-\sigma}\right)$$ se le conoce como [**logit**](https://en.wikipedia.org/wiki/Logit) y modela el log momio: $$\log \left(\frac{p(\mathcal{C}_1|x)}{p(\mathcal{C}_0|x)}\right)=\log \left(\frac{p(\mathcal{C}_1|x)}{1-p(\mathcal{C}_1|x)}\right)$$que tiene una interpretación directa en términos de las probabilidades de pertenencia a cada clase $\mathcal{C}_0,\mathcal{C}_1$.  ¿Modelo en regresión logística de dos clases? De forma similar como en el modelo por mínimos cuadrados lineales (ver [4.3.Minimos_cuadrados_R](https://github.com/ITAM-DS/analisis-numerico-computo-cientifico/blob/master/temas/IV.optimizacion_convexa_y_machine_learning/4.3.Minimos_cuadrados_R.ipynb)) 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: $$\beta^T x+\beta_0=a(x|\beta_0,\beta)=\log \left(\frac{p(\mathcal{C}_1|x)}{p(\mathcal{C}_0|x)}\right).$$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: $$p(\mathcal{C}_1 | x ) = y(x | \beta_0, \beta) = \frac{1}{1+ e^{-(\beta_0, \beta)^T x}}$$que se lee "la probabilidad de pertenencia a la clase $\mathcal{C}_1$ dado el vector de atributos $x$ es igual a $y$". **Comentarios:** * El modelo con $2$ parámetros $\beta_0, \beta_1$ se ve como:$$p(\mathcal{C}_1 | x ) = y(x | \beta_0, \beta) = \frac{1}{1+ e^{-(\beta_0 + \beta_1x)}}$$con $x \in \mathbb{R}$. * El modelo puede extenderse utilizando $n+1$ funciones conocidas $\phi_j:\mathbb{R} \rightarrow \mathbb{R}$, $\phi_j(x)$ $j=0,\dots, n$ por lo que si $\phi(x)=(\phi_0(x),\phi_1(x),\dots,\phi_n(x))^T$ y $\beta_0 \in \mathbb{R}$, $\beta \in \mathbb{R}^n$, entonces se tiene el modelo por regresión logística:$$p(\mathcal{C}_1|\phi(x))=y(x|\beta_0, \beta)= \frac{1}{1+ e^{-(\beta_0, \beta)^T \phi(x)}}$$ * 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}$.  ¿Cómo se ajustan los parámetros del modelo por regresión logística de dos clases? 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](https://en.wikipedia.org/wiki/Likelihood_function):$$\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: $$\displaystyle \max_{(\beta_0, \beta) \in \mathbb{R}^{n+1}} \mathcal{L}(\beta_0, \beta|x)=\displaystyle \prod_{i=0}^n y_i^{\hat{y}_i}(1-y_i)^{1-\hat{y}_i}$$ 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**](https://en.wikipedia.org/wiki/Deviance_(statistics): $$\begin{eqnarray}\displaystyle \min_{(\beta_0, \beta) \in \mathbb{R}^{n+1}}\mathcal{D}(\beta_0, \beta|x)&=&-2\ell(\beta_0, \beta|x) \nonumber \\&=&2\displaystyle \sum_{i=1}^m\log(1+\exp((\beta_0, \beta)^Tx_i))-\hat{y}_i(\beta_0, \beta)^Tx_i \nonumber\end{eqnarray}$$ **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(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: $$A = \left[\begin{array}{c}x_0\\x_1\\\vdots\\x_m\end{array}\right]=\left[\begin{array}{cccc}x_{01} & x_{02}&\dots& x_{0n}\\x_{11}& x_{12}&\dots& x_{1n}\\\vdots &\vdots& \vdots&\vdots\\x_{n1} &x_{n2}&\dots&x_{nn}\\\vdots &\vdots& \vdots&\vdots\\x_{m1} &x_{m2}&\dots&x_{mn}\end{array}\right] \in \mathbb{R}^{(m+1)x(n+1)}$$ \*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: $$p=\left[\begin{array}{c}y_0(x_0|\beta_0,\beta)\\y_1(x_1|\beta_0,\beta)\\\vdots \\y_m(x_m|\beta_0,\beta)\end{array}\right]=\left[\begin{array}{c}p(\mathcal{C}_1|x_0)\\p(\mathcal{C}_1|x_1)\\\vdots \\p(\mathcal{C}_1|x_m)\end{array}\right] \in \mathbb{R}^{n+1},\hat{y}=\left[\begin{array}{c}\hat{y}_0\\\hat{y}_1\\\vdots \\\hat{y}_m\end{array}\right] \in \mathbb{R}^{m+1}$$ 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](https://github.com/ITAM-DS/analisis-numerico-computo-cientifico/blob/master/temas/IV.optimizacion_convexa_y_machine_learning/4.1.Optimizacion_numerica_y_machine_learning.ipynb). 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](https://github.com/ITAM-DS/analisis-numerico-computo-cientifico/blob/master/temas/IV.optimizacion_convexa_y_machine_learning/4.2.Algoritmos_para_UCO.ipynb).  Ejemplos Como primer ejemplo utilizamos el conocido *dataset* de [*Iris*](https://en.wikipedia.org/wiki/Iris_flower_data_set) en el que se muestran tres especies del género *Iris*. Las especies son: [*I. setosa*](https://en.wikipedia.org/wiki/Iris_setosa), [*I. virginica*](https://en.wikipedia.org/wiki/Iris_virginica) y [*I. versicolor*](https://en.wikipedia.org/wiki/Iris_versicolor):Imagen obtenida de [Iris Dataset](https://rpubs.com/AjinkyaUC/Iris_DataSet).<jupyter_code>library(datasets)
library(dplyr)
library(glmnet)
data(iris)
summary(iris)

print("Número de renglones y columnas:")
nrow(iris)
ncol(iris)<jupyter_output>[1] "Número de renglones y columnas:"<jupyter_text>**Para el modelo de regresión logística en dos clases utilizamos a la clase virginica y versicolor.**<jupyter_code>iris_subset <- filter(iris, Species %in% c("virginica", "versicolor"))

head(iris_subset)
nrow(iris_subset)<jupyter_output><empty_output><jupyter_text>**Obsérvese que la variable `Sepal.Length` y `Petal.Length` están altamente correlacionadas por lo que quitaremos a `Petal.length`:**<jupyter_code>cor(as.matrix(iris_subset[,-ncol(iris_subset)]))
A<-iris_subset[,-c(3,ncol(iris_subset))]
head(A)<jupyter_output><empty_output><jupyter_text>La variable $\hat{y}$ estará dada por las clases `versicolor` y `virginica`:<jupyter_code>head(iris_subset[,ncol(iris_subset)])
tail(iris_subset[,ncol(iris_subset)])<jupyter_output><empty_output><jupyter_text>**Transformamos a tipo numérica la columna `Species`:**<jupyter_code>iris_subset$Clase <- as.numeric(iris_subset$Species %in% 'versicolor')
head(iris_subset)
tail(iris_subset)<jupyter_output><empty_output><jupyter_text>La clase $\mathcal{C}_1$ es `versicolor` y $\mathcal{C}_0$ es `virginica`. **Utilizamos el paquete la función [glm](https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/glm) del paquete  [stats](https://www.rdocumentation.org/packages/stats/versions/3.6.2) para realizar el ajuste:**<jupyter_code>glm.out <- glm(Clase ~ Sepal.Length + Sepal.Width  + Petal.Width,
               data = iris_subset,family = binomial)

summary(glm.out)<jupyter_output><empty_output><jupyter_text>Y nos interesan las estimaciones de los coeficientes $\beta_0, \beta_1, \beta_2, \beta_3$:<jupyter_code>print(glm.out$coefficients)<jupyter_output>(Intercept) Sepal.Length  Sepal.Width  Petal.Width 
   20.287329    -1.295129     4.823274   -15.922660<jupyter_text>Estimación vía método descenso con la dirección de Newton La función objetivo como se revisó en la sección anterior está dada por la expresión de la devianza: $$\mathcal{D}(\beta_0, \beta|x)=-2\ell(\beta_0, \beta|x) = 2\displaystyle \sum_{i=1}^m\log(1+\exp((\beta_0, \beta)^Tx_i))-\hat{y}_i(\beta_0, \beta)^Tx_i$$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$:**<jupyter_code>A<-cbind(rep(1,nrow(A)),A)
head(A)
y_hat <- iris_subset$Clase<jupyter_output><empty_output><jupyter_text>La clase $\mathcal{C}_1$ es `versicolor`:<jupyter_code>head(y_hat)<jupyter_output><empty_output><jupyter_text>La clase $\mathcal{C}_0$ es `virginica`:<jupyter_code>tail(y_hat)<jupyter_output><empty_output><jupyter_text>**Cargamos las funciones programadas en el directorio `algoritmos/R`:**<jupyter_code>#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=""))<jupyter_output><empty_output><jupyter_text>**Convertimos a tipo matriz:**<jupyter_code>A<-as.matrix(A)
head(A)<jupyter_output><empty_output><jupyter_text>**Función objetivo:** $$2\displaystyle \sum_{i=1}^m\log(1+\exp((\beta_0, \beta)^Tx_i))-\hat{y}_i(\beta_0, \beta)^Tx_i$$<jupyter_code>fo <- function(beta){
    mat_vec <- A%*%beta
    2*sum(log(1+exp(mat_vec))-y_hat*mat_vec)
    }<jupyter_output><empty_output><jupyter_text>**Punto inicial $\beta ^{(0)}:$**<jupyter_code>beta_0 <- c(1,1,1,1)<jupyter_output><empty_output><jupyter_text>**$\beta^*$:**<jupyter_code>beta_ast_glm <- glm.out$coefficients

print(beta_ast_glm)
p_ast <- fo(beta_ast_glm)
p_ast<jupyter_output><empty_output><jupyter_text>**Argumentos para el método de Newton:**<jupyter_code>tol <- 1e-8
tol_backtracking <- 1e-14
maxiter <- 30
l<-Newtons_method(fo, beta_0, tol, tol_backtracking, beta_ast_glm, p_ast, maxiter)
beta <- l[[1]]
total_of_iterations <- l[[2]]
Err_plot <- l[[3]]
beta_plot <- l[[4]]<jupyter_output><empty_output><jupyter_text>**Secuencia de minimización $\beta^{(k)}$:**<jupyter_code>beta_plot<jupyter_output><empty_output><jupyter_text>**$\beta$ aproximada por el método de Newton:**<jupyter_code>print(beta)<jupyter_output>[1]  20.262418  -1.292617   4.822189 -15.915167<jupyter_text>**$\beta^*$:**<jupyter_code>print(beta_ast_glm)
compute_error(x_obj = beta_ast_glm,x_approx = beta)<jupyter_output><empty_output><jupyter_text>**Tenemos alrededor de 3 dígitos de precisión.** **Comparemos con la función [glmnet](https://www.rdocumentation.org/packages/glmnet/versions/3.0-2/topics/glmnet) del paquete del mismo nombre (no usamos intercepto):**<jupyter_code>fit<-glmnet(A[,-1],y_hat,family="binomial",alpha=0,lambda=0,standardize=F,nlambda=1, thresh=1e-8)

beta_ast_glmnet <- as.matrix(fit$beta)
print(beta_ast_glmnet)<jupyter_output>s0
Sepal.Length  -1.295005
Sepal.Width    4.823018
Petal.Width  -15.922371<jupyter_text>**Error relativo:**<jupyter_code>compute_error(x_obj = beta_ast_glmnet,x_approx = beta[-1])<jupyter_output><empty_output><jupyter_text>**Tenemos alrededor de 3 dígitos de precisión.**  Cálculo de probabilidades de pertenencia a las clases $\mathcal{C}_0 =$ `virginica`, $\mathcal{C}_1 =$ `versicolor` **Para individuo $i$ se tiene:** $$p(\mathcal{C}_1|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.$$ **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`**:<jupyter_code>1/(1+exp(-sum(A[1,]*beta_ast_glm)))<jupyter_output><empty_output><jupyter_text>**Estimación de la probabilidad de pertenencia a la clase $\mathcal{C}_0$: `virginica` con lo estimado por `glm`**:<jupyter_code>exp(-sum(A[1,]*beta_ast_glm))/(1+exp(-sum(A[1,]*beta_ast_glm)))<jupyter_output><empty_output><jupyter_text>**Estimación de la probabilidad de pertenencia a la clase $\mathcal{C}_1$:`versicolor` con lo estimado por descenso por Newton**:<jupyter_code>1/(1+exp(-sum(A[1,]*beta)))<jupyter_output><empty_output><jupyter_text>**Estimación de la probabilidad de pertenencia a la clase $\mathcal{C}_0$:`virginica` con lo estimado por descenso por Newton**:<jupyter_code>exp(-sum(A[1,]*beta))/(1+exp(-sum(A[1,]*beta)))<jupyter_output><empty_output><jupyter_text>**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`**:<jupyter_code>1/(1+exp(-sum(A[nrow(A),]*beta_ast_glm)))<jupyter_output><empty_output><jupyter_text>**Estimación de la probabilidad de pertenencia a la clase $\mathcal{C}_0$: `virginica` con lo estimado por `glmn`**:<jupyter_code>exp(-sum(A[nrow(A),]*beta_ast_glm))/(1+exp(-sum(A[nrow(A),]*beta_ast_glm)))<jupyter_output><empty_output><jupyter_text>**Estimación de la probabilidad de pertenencia a la clase $\mathcal{C}_1$:`versicolor` con lo estimado por descenso por Newton**:<jupyter_code>1/(1+exp(-sum(A[nrow(A),]*beta)))<jupyter_output><empty_output><jupyter_text>**Estimación de la probabilidad de pertenencia a la clase $\mathcal{C}_0$:`virginica` con lo estimado por descenso por Newton**:<jupyter_code>exp(-sum(A[nrow(A),]*beta))/(1+exp(-sum(A[nrow(A),]*beta)))<jupyter_output><empty_output><jupyter_text>**A continuación consideramos regularización con ridge y lasso, ver [4.3.Minimos_cuadrados_R.ipynb](https://github.com/ITAM-DS/analisis-numerico-computo-cientifico/blob/master/temas/IV.optimizacion_convexa_y_machine_learning/4.3.Minimos_cuadrados_R.ipynb).**  Regularización ridge con intercepto<jupyter_code>reg<-0.5
fit<-glmnet(A,y_hat,family="binomial",alpha=0,lambda=reg,standardize=F,nlambda=1, thresh=1e-8)<jupyter_output><empty_output><jupyter_text>**$\beta^*:$**<jupyter_code>beta_ast<-as.matrix(fit$beta)
print(beta_ast)<jupyter_output>s0
rep(1, nrow(A))  0.00000000
Sepal.Length    -0.24479206
Sepal.Width     -0.07182931
Petal.Width     -0.30070071<jupyter_text>**Función objetivo de la función glmnet con regularización *ridge*:** $$\frac{1}{m}\displaystyle \sum_{i=1}^m\log(1+\exp((\beta_0, \beta)^Tx_i))-\hat{y}_i(\beta_0, \beta)^Tx_i + \frac{\lambda}{2}||\beta||_2^2$$<jupyter_code>mpoints<-nrow(A)
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)])
    }<jupyter_output><empty_output><jupyter_text>**Punto inicial $\beta ^{(0)}:$**<jupyter_code>beta_0 <- c(0,0,0,0)
p_ast <- fo(beta_ast)
p_ast
tol <- 1e-12
tol_backtracking <- 1e-14
maxiter <- 30

l<-Newtons_method(fo, beta_0, tol, tol_backtracking, beta_ast, p_ast, maxiter)
beta <- l[[1]]
total_of_iterations <- l[[2]]
Err_plot <- l[[3]]
beta_plot <- l[[4]]<jupyter_output><empty_output><jupyter_text>**$\beta$ aproximada por el método de Newton:**<jupyter_code>print(beta)
print(beta_ast)<jupyter_output>s0
rep(1, nrow(A))  0.00000000
Sepal.Length    -0.24479206
Sepal.Width     -0.07182931
Petal.Width     -0.30070071<jupyter_text>**Error relativo:**<jupyter_code>compute_error(beta_ast[-1],beta[-1])<jupyter_output><empty_output><jupyter_text>**Sin considerar el intercepto tenemos alrededor de $6$ dígitos de precisión.**  Regularización Ridge sin intercepto<jupyter_code>A<-A[,-1]
head(A)<jupyter_output><empty_output><jupyter_text>**Solución con glmnet:**<jupyter_code>fit<-glmnet(A,y_hat,family="binomial",alpha=0,lambda=reg,standardize=F,nlambda=1, intercept=F, thresh=1e-8)

beta_ast <- as.matrix(fit$beta)<jupyter_output><empty_output><jupyter_text>**$\beta^*:$**<jupyter_code>print(beta_ast)<jupyter_output>s0
Sepal.Length  0.02157718
Sepal.Width   0.05710003
Petal.Width  -0.24262286<jupyter_text>**Solución vía método de Newton:** **Punto inicial $\beta ^{(0)}:$**<jupyter_code>beta_0 <- c(0,0,0)
p_ast <- fo(beta_ast)
p_ast
tol <- 1e-12
tol_backtracking <- 1e-14
maxiter <- 30
fo <-function(beta){
    mat_vec <- A %*% beta
    1/mpoints*(sum(log(1 + exp(mat_vec)) - y_hat * mat_vec)) + reg/2*sum(beta*beta)
    }
l<-Newtons_method(fo, beta_0, tol, tol_backtracking, beta_ast, p_ast, maxiter)
beta <- l[[1]]
total_of_iterations <- l[[2]]
Err_plot <- l[[3]]
beta_plot <- l[[4]]<jupyter_output><empty_output><jupyter_text>**$\beta$ aproximada por el método de Newton:**<jupyter_code>print(beta)
print(beta_ast)<jupyter_output>s0
Sepal.Length  0.02157718
Sepal.Width   0.05710003
Petal.Width  -0.24262286<jupyter_text>**Error relativo:**<jupyter_code>compute_error(x_approx = beta, x_obj = beta_ast)<jupyter_output><empty_output><jupyter_text>**Tenemos alrededor de $5$ dígitos de precisión.**  Regularización lasso sin intercepto **Función objetivo de la función glmnet con regularización *lasso*:** $$\frac{1}{m}\displaystyle \sum_{i=1}^m\log(1+\exp((\beta_0, \beta)^Tx_i))-\hat{y}_i(\beta_0, \beta)^Tx_i + \lambda \alpha||\beta||_1$$<jupyter_code>reg <- 0.01<jupyter_output><empty_output><jupyter_text>**Solución vía `glmnet`:**<jupyter_code>fit<-glmnet(A,y_hat,family="binomial",alpha=1,lambda=reg,standardize=F,nlambda=1, intercept=F, thresh=1e-8)

beta_ast <- as.matrix(fit$beta)<jupyter_output><empty_output><jupyter_text>**$\beta^*:$**<jupyter_code>print(beta_ast)<jupyter_output>s0
Sepal.Length  0.6246049
Sepal.Width   3.0968033
Petal.Width  -7.8624767<jupyter_text>**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:$$\text{quita_signo}(x) = \begin{cases}x & \text{si } x > 0,\\-x & \text{si } x <0,\\\approx 2.22 \times 10^{-308} & \text{si } |x| \approx 0\end{cases}$$ **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](https://github.com/ITAM-DS/analisis-numerico-computo-cientifico/blob/master/temas/I.computo_cientifico/1.2.Sistema_de_punto_flotante.ipynb).**<jupyter_code>quita_signo<-function(beta){
    beta<-sign(beta)*beta
    ind <- beta < .Machine$double.xmin & beta > -.Machine$double.xmin 
    beta[ind] <- .Machine$double.xmin 
    beta
    }<jupyter_output><empty_output><jupyter_text>**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](https://en.wikipedia.org/wiki/Subderivative), útiles para generalizar las derivadas para funciones convexas que no son diferenciables en puntos de su dominio.<jupyter_code>fo <-function(beta){
    mat_vec <- A %*% beta
    1/mpoints*(sum(log(1 + exp(mat_vec)) - y_hat * mat_vec)) + reg*sum(quita_signo(beta))
    }<jupyter_output><empty_output><jupyter_text>**Punto inicial $\beta ^{(0)}:$**<jupyter_code>beta_0<-c(1,1,1)
l<-Newtons_method(fo, beta_0, tol, tol_backtracking, beta_ast, p_ast, maxiter)
beta <- l[[1]]
total_of_iterations <- l[[2]]
Err_plot <- l[[3]]
beta_plot <- l[[4]]<jupyter_output><empty_output><jupyter_text>**$\beta$ aproximada por el método de Newton:**<jupyter_code>print(beta)
print(beta_ast)<jupyter_output>s0
Sepal.Length  0.6246049
Sepal.Width   3.0968033
Petal.Width  -7.8624767<jupyter_text>**Error relativo:**<jupyter_code>compute_error(x_approx = beta, x_obj = beta_ast)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17580680732191428, 'fraction_numerical': 0.02870690288906859, 'mean_word_length': 5.794955044955045, 'pattern_counts': {'":': 0, '<': 282, '<?xml version=': 0, '>': 206, 'https://': 30, 'lorem ipsum': 0, 'www.': 5, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 2, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '53', 'dir': 'jupyter-structured-clean-dedup', 'id': '383573', 'license': 'Apache-2.0', 'n_tokens_mistral': 11294, 'n_tokens_neox': 10767, 'n_words': 2089, 'path': 'temas/IV.optimizacion_convexa_y_machine_learning/4.3.Regresion_logistica_R.ipynb', 'repo_name': '123972/analisis-numerico-computo-cientifico', 'size': '27202'}