Algoritmos de Machine Learning en R project

“The best thing about R is  that it was developed by statisticians. The worst thing about R is that….it was developed by statisticians.”  Bo Cowgill, Google, Inc.

Todas las entradas de esta categoría muestran algoritmos y técnicas clasificadas como algoritmos de Machine Learning o Aprendizaje Estadístico. En cada entrada comparto el código  y algunos repositorios donde pueden bajar datos y notas sobre los temas.

Las ventajas de usar R project para dar estos ejemplos son las siguientes:

1.-R project es software libre u open source, por lo cual cualquier persona pude descargarlo.

2.-El código que comparto solo debe de ser copiado y pegado en un Script o en la consola  de R project para reproducir los resultados.

Otra ventaja de aprender estas técnicas con ejemplos en R, es que muchos de los algoritmos son acompañados de alguna técnica gráfica para visualizar el resultado. R project por default cuenta con una gran cantidad de herramientas gráficas o se puede hacer uso de librerías como ggplot2, lattice, ggvis para hacer gráficas aún más elaboradas y visualmente más agradables.

El entendimiento de la técnica después de realizar un ejemplo desde mi perspectiva genera interés por saber como aplicarla a otros datos y nos invita a conocer más sobre el algoritmo desde un punto de vista no solo práctico, sino también teórico.

Evité en la medida de lo posible hacer uso de ecuaciones y en cada ocasión que se menciona algún concepto de matemáticas trato de explicarlo de un modo sencillo. Si este concepto requiere a mi parecer una mejor lectura, agregue algunas ligas para consultar en la red.

La única meta que me propuse fue explicar las cosas de modo tal que cualquier persona con un poco de curiosidad pueda entender algo de cada entrada.

Teniendo en mente lo anterior, espero que  el objetivo de trasmitir el conocimiento desde un ejemplo concreto haya sido logrado. En caso de que no sea así, sería grato saber qué parte no es clara y como mejorar este material.

Resumen de las entradas.

El orden siguiente es el que considero apropiado para revisar cada entrada.

1.-¿Qué se necesita conocer de R project para aprender algunas técnicas de Machine Learning?

Es una breve explicación sobre R project y sobre algunas de las funciones con las que cuenta por default. También menciono las librería que se requieren para replicar los ejemplos, aunque en cada entrada menciono algo breve sobre la librerías requeridas.

2.-Análisis Confirmatorio vs Exploratorio

Explico la diferencia entre los dos tipos de análisis y cómo están relacionados. Clasifico el tipo de medidas que uno puede analizar en los datos y pongo algunos ejemplos del tipo de gráficas que pueden ayudarnos a visualizar las medidas al momento de explorar la información.

3.-Regresión Lineal, un ejemplo sencillo.

Esta técnica es muy sencilla de implementar y los datos que nos regresa el comando en R debe ser entendidos para poder evaluar la calidad del resultado obtenido. Omito detalles formales del tema y comparto algunos ejemplos de regresiones simples y múltiples. Es un tema con muchos detalles importantes, por ejemplo la revisión de las posibles fallas al construir un modelo lineal o multilineal. No toco mucho  esos temas, pero en las referencias pueden encontrar observaciones y técnicas para conocer más al respecto.

4.-Regresión no lineal, Cross-Validation y Regularization. 

Los tres temas pueden haber requerido una entrada para cada uno, pero funciona bien mostrar como se relacionan. El objetivo es mostrar la implementación de la regresión no lineal y validar que el modelo elegido es el adecuado, para evitar la sobre estimación de los datos. En esto último las cross-validation y la regularización ayudan a tener un método para hacer una elección del modelo. Se dan algunos ejemplos del uso de los tres conceptos y técnicas.

5.-Clasificación binaria….Naive Bayes. 

Es uno de los algoritmos más usados por sencillo y por que puede ser modificado para tener una clasificación rápida en varias categorías. Doy el ejemplo clásico sobre detección de Spam y trato de complementarlo con el uso de una librería en la cual se cuenta con una función para estimar el modelo. Agregué una pequeña explicación sobre la probabilidad condicional y sobre el uso de la librería tm para analizar textos en R project. Espero ser claro y que permita comprenderse el ejemplo.

6.-Análisis de Cluster, un ejemplo sencillo. 

Se explica de modo breve como funciona la detección de Clusters. Como el tema requiere hacer uso de la noción de distancia o métrica, trato de explicar el concepto. Omití mucho material, sobre todo en el uso de distintas métricas para distintos tipos de variables.

7.-Un ejemplo sencillo sobre Análisis de Componentes Principales (PCA, Principal Component Analysis).

La técnica es muy usada y conocida para hacer reducción de dimensiones, explico qué es reducir dimensiones y como interpretar los resultado de la técnica. Como ejemplo principal hago una reducción de 25 indicadores para formar uno que resuma la información del resto. Replico un ejemplo de análisis de la digitalización de los números escritos a mano, este es un ejemplo clásico y la base de datos fue objeto de estudio para comparar técnicas de clasificación.

8.-Sobre sistemas de recomendación, ejemplo sencillo. 

Hago un breve ejemplos de un sistema de recomendación usando la técnica de k-vecinos cercanos (K-nn), esta técnica es no-paramétrica. Para ejemplificar como funciona pongo un ejemplo inicial y trato de explicar qué realiza el algoritmo y aplicarlo a un caso concreto. Los sistemas de recomendación prácticamente dan la libertad de elegir la técnica de clasificación que uno considere adecuada para el tipo de datos o servicio, por lo cual existen libertad de replicar este ejemplo haciendo uso de otro técnicas de clasificación.

9.-Optimización.

Esta entrada creo que es extensa y puede ser un poco amplia,ya que traté de explicar desde aspectos básicos de mínimo y máximos hasta algoritmo estocásticos y genéticos. Puse un ejemplo de como usar el código y en los casos más sofisticados, la parte estocástica, comparto un ejemplo que permitiera comparar los resultado obtenidos con diversos tipos de algoritmos, no puse los tiempos que tardan en dar la solución pero en caso de replicarlos será notorio el costo del procesamiento.

10.-Máquina de soporte Vectorial (SVM-Sopport Vector Machine).

La técnica sumamente usada y funciona bien para datos no-lineales, se dan varios ejemplos de como usar dos librerías que permiten calcular SVM con distintos kernels. En uno de los ejemplos de hace la comparación de varios modelos para determinar por medio de MSE cual es el mejor. Son ejemplos muy sencillos y trato de agregar explicaciones a cada ejemplo.

Lo que no está aún en estas entradas, pero que estará en inicio de Marzo 2016.

Dejé fuera por el momento varios temas fuera de las entradas, ejemplo Redes neuronales , Arboles de decisión, Random Forests, AdaBoost. Espero posteriormente agregar entradas sobre esos algoritmos.

En general todos las entradas tratan de simplificar y resumir lo que se puede encontrar en el texto  “Machine Learning for Hackers“, libro escrito por Jhon M. White y Drew Conway. Si bien, el libro en si es sencillo y requiere un bajo nivel de matemáticas para estudiarse, traté de hacer más sencillas las explicaciones y de complementarlo con otros ejemplos.

La recomendación es complementar las explicaciones de las entradas con la lectura del texto.  Jhon M White compartió su código  y datos usados en los ejemplos en el repositorio Github y pueden descargar el código presionando estas letras en rojo.

Algunos de los códigos deben de revisarse ya que las librerías que usaron en el año de la edición del texto ya fueron actualizadas, por lo cual es generan algunos errores.

Dos temas que vienen en el texto y decidí omitir son generación de un rango y análisis de redes sociales o de grafos. En otro momento comparto ejemplos respecto al tema.

El segundo es sumamente interesante; análisis de redes sociales, pero el código con el que cuenta el texto ya no es funcional. Así que se requiere rehacer buena parte del código y diseñar nuevos ejemplos, para eso se necesita cierto tiempo. Más aún debido al boom de las redes sociales existen otro tipo de herramientas sobre las cuales se puede escribir.

Para tener una idea de lo nuevo e interesante se pueden consultar las siguientes ligas para echar un vistazo:

Deseo que el material que se encuentra en esta categoría les sea útil y sobre todo les deje con ganas de aprender más al respecto y no solo desde el lado aplicado, sino también desde el lado teórico.

En la siguiente cita aplique “instrumento= algoritmo”.

“Tienes que aprender a tocar tu instrumento. Después debes practicar, practicar y practicar. Y después, cuando finalmente estás en el escenario, olvídalo todo y ulula.” Charlie Parker

 Nota: Todas las críticas y comentarios respecto a las entradas son bienvenidos.
Anuncios

Máquina de soporte Vectorial (SVM-Sopport Vector Machine)

En la entrada explico la técnica de máquina de soporte vectorial (SVM-Sopport Vector Machine),  es un conjunto de algoritmos que pueden aplicarse a problemas de  clasificación o regresión.

La explicación no es exhaustiva ni mucho menos menciono aspectos teóricos de la técnica, tan solo trato de mostrar con varios ejemplos como funciona.

Existe mucho material para revisar los detalles y como menciona el autor de la referencia [2], existen ahora nuevas técnicas para clasificar pero es casi obligado en Machine Learning conocer el algoritmo SVM por la relevancia que han tenido.

Una de esas nuevas  técnicas es Relevance Vector Machine (RVM), la cual es una variación de SVM en su versión Bayesiana y que permite ver la técnicas SVM desde un punto de vista probabilista.

En otras dos entradas comenté sobre los temas  Naive Bayes y Regresión Lineal, la primera técnica es por excelencia una técnica de clasificación y la segunda una de predicción pero se puede desarrollar o modificar para usarla para clasificar.

Las técnicas de SVM cubren ambos aspectos, clasificación y predicción. En la entrada solo comento sobre como usar esta técnica para hacer clasificaciones y muestro un ejemplo sencillo de como usarla para hacer predicción, lo malo de esta familia de técnica es que al construir un modelo uno desea o busca en algunas ocasiones tener una ecuación explicita para analizar el comportamiento de las variables. En este caso eso no es posible, pero para ciertas situaciones no es necesario saber como afectan las variaciones de las variables predictoras, ejemplo para hacer algo de Marketing uno quizás solo busca la clasificación mejor de la muestra o para clasificar textos, quizás no interesa la forma del modelo porque no tenemos una variable respuesta a modelar sino solo entras para clasificar.

Para los ejemplos hago uso de dos librerías en R que calculan los algoritmos de SVM; e1071 y kernlab. También se puede usar para realizar SVM como método de clasificación para dos clases la librería svmpath.

Algunos ejemplos con las librerías

Antes de hacer el ejemplo del texto Machine Learning for Hackers[1] , hago uso de la librería e1071 para mostrar uno de los conceptos claves en SVM, lo que es un hiper-plano. En este caso lo “hiper” se reduce a una recta en el plano. El código es el siguiente:

#Código de muestra de la técnica SVM
#Se carga a librería e1071
library(e1071)
set.seed(1)
#Se generan los datos
x=matrix(rnorm(100*2), ncol=2)
y=c(rep(-1,10), rep(1,10))
x[y==1,]=x[y==1,] + 1
#Se hace la gráfica de los puntos generados en el plano
plot(x, col=(3-y),main="Datos mezclados",xlab="Eje X",ylab="Eje y")
#Se construye un Data.Frame
dat=data.frame(x=x, y=as.factor(y))

#Se calcula SVM para el kernel lineal
svmfit=svm(y~., data=dat, kernel="linear", cost=10)
plot(svmfit, dat)
#Se puede revisar los valores y el resumen de los estimado en R project para el modelo SVM
#svmfit$index
#summary(svmfit)

#Se estima el modelo SVM con kernel polinomial
svmfit=svm(y~., data=dat, kernel="polynomial", cost=10)
plot(svmfit, dat)
#Se  estima el modelo SVM con el kernel sigmoid
svmfit=svm(y~., data=dat, kernel="sigmoid", cost=10)
plot(svmfit, dat)
#Se estima el modelo SVM con el kernel radial
svmfit=svm(y~., data=dat, kernel="radial", cost=10)
plot(svmfit, dat)

Cuando se corre el código se obtienen cuatro gráficas, la primera solo presenta los datos a clasificar en el plano, la segunda, tercera, cuarta y quinta gráfica muestran la clasificación en dos clases por medio de variaciones al modelo SVM.

Datos_meclados_en_el_planoSVM_Lineal

SVM_polinomial

SVM_sigmoid

SVM_Radial

Cada gráfica muestra una separación del conjunto de puntos en dos grupos, pese a que no tienen un título cada gráfica que indique el kernel que se usó, la secuencia es correcta con respecto al código. Entonces lo que hace el tipo de kernel es definir un tipo de línea o frontera entre las clases. Se observa que cuando el kernel es lineal se obtienen como frontera entre las clases como la recta de un ajuste lineal. Esa rectas, para el caso lineal la recta es lo que se llama hiperplano.

Para comparar el funcionamiento del algoritmo como técnica de clasificación, hago una comparación con el método Naive Bayes aplicándolo a una muestra 4601 correos electrónico analizados y con una clasificación de ser spam o no. Estos datos se encuentran en la librería kernlab y hago uso de la función naivebayes de la librería e1071.

Los datos son 4601 correos con 58 variables asignadas, es decir; la tabla de datos tienen 4601 renglones y 58 columnas. La última columna indica si el correo fue clasificado como spam o no, esto permite comparar al final la clasificación de una muestra. En este caso se toman 100 correos para evaluar el método SVM kernel radial y polinomial.

#Comparación con otro método de clasificación
library(kernlab)
data(spam)
index=sample(1:nrow(spam),100)
fm=ksvm(type~.,data=spam[-index,],kernel="rbfdot",kpar="autolibrary(e1071)matic",C=60,cross=5)
pred=predict(fm,spam[index,])
table(pred,spam[index,58])
#Los datos que se obtienen al ejecutar la clasificación son:
#pred nonspam spam
#nonspam  62   2
#spam      2   34

fm=ksvm(type~.,data=spam[-index,],kernel="polydot",kpar="automatic",C=60,cross=5)
pred=predict(fm,spam[index,])
table(pred,spam[index,58])
#Los datos que se obtienen al ejecutar la clasificación son:
#pred nonspam spam
#nonspam 62    4
#spam     2    32
#Comparación con otro método de clasificación
library(e1071)
model=naiveBayes(type~.,data=spam[-index,])
pred=predict(mol,spam[index,])
table(pred,spam[index,58])
#pred nonspam spam
#nonspam 38    1
#spam    26    35

Calculando el porcentaje de eficiencia de los métodos, SVM tienen 96%  y 94% de eficiencia, por otro lado el método Naive Bayes tienen 73% de aciertos. En este ejemplo se observa que el método es eficiente para hacer clasificaciones, en este último ejemplo se tienen 58 variables no todas siguen comportamientos lineales, por lo cual tomar un kernel no lineal para  SVM  resulta apropiado.

Lo que se observa es que SVM puede servir como un método de clasificación para dos clases, pero también es adecuado para construir un modelo para clasificar para más de dos clases.

El  ejemplos del libro Machine Learning for Hackers[1], del cual se pueden descargar los datos en github, inicialmente muestra un comparativo entre un el método de clasificación lineal logit contra SVM.

Hago la regresión logística para los datos que cuentan con una variable indicadora para distinguir entre aceptados y no aceptados y se compara la estimación con el método SVM.

#Cargamos los datos
df <- read.csv(file.path('data', 'df.csv'))

#Calculamos la regresión logística
logit.fit <- glm(Label ~ X + Y,
 family = binomial(link = 'logit'),
 data = df)

#Estimamos la predicción de los aceptados
logit.predictions <- ifelse(predict(logit.fit) > 0, 1, 0)

#Estimamos la media de la predicción

mean(with(df, logit.predictions == Label))
#[1] 0.5156

 El cálculo de SVM se hace mediante la librería e1071.

#Cargamos la librería

library('e1071')

#Estimamos el modelo por medio de SVM

svm.fit <- svm(Label ~ X + Y, data = df)
svm.predictions <- ifelse(predict(svm.fit) > 0, 1, 0)

#Calculamos la media de las predicciones
mean(with(df, svm.predictions == Label))
#[1] 0.7204

 Los resultados son los siguientes, el método lineal estima correctamente el 51%  de los datos y  SVM el 72% .

Para este ejemplo el kernel que se usa para SVM es lineal  y sin embargo arroja mejores resultados que la regresión logística.

La gráfica de los resultados se verían así:

#Graficamos los resultados
library("reshape")
df <- cbind(df,
 data.frame(Logit = ifelse(predict(logit.fit) > 0, 1, 0),
 SVM = ifelse(predict(svm.fit) > 0, 1, 0)))

predictions <- melt(df, id.vars = c('X', 'Y'))

ggplot(predictions, aes(x = X, y = Y, color = factor(value))) +
 geom_point() +
 facet_grid(variable ~ .)+ggtitle("Ejemplo SVM vs Logit") + theme(plot.title = element_text(lineheight=.9, face="bold"))

Ejemplo_1

En las primeras cuatro gráficas de la entrada mostré que cuando se cambia el tipo de kernel para método SVM las líneas de las clases mediante las cuales separa los datos se comportan de modo distinto. Al aplicar varios tipos de kernel para los datos df, se muestra un gráfico comparativo entre el comportamiento de las fronteras entre los datos.

#Comparamos distintos kernel
#Usamos los datos anteriores

df <- df[, c('X', 'Y', 'Label')]

#Calculamos SVM con el Kernel Lineal
linear.svm.fit <- svm(Label ~ X + Y, data = df, kernel = 'linear')
with(df, mean(Label == ifelse(predict(linear.svm.fit) > 0, 1, 0)))

#Calculamos SVM con el Kernel Polinimial
polynomial.svm.fit <- svm(Label ~ X + Y, data = df, kernel = 'polynomial')
with(df, mean(Label == ifelse(predict(polynomial.svm.fit) > 0, 1, 0)))

#Calculamos SVM con el Kernel Radial o RBF 
radial.svm.fit <- svm(Label ~ X + Y, data = df, kernel = 'radial')
with(df, mean(Label == ifelse(predict(radial.svm.fit) > 0, 1, 0)))

#Calculamos SVM con el Kernel Sigmoid
sigmoid.svm.fit <- svm(Label ~ X + Y, data = df, kernel = 'sigmoid')
with(df, mean(Label == ifelse(predict(sigmoid.svm.fit) > 0, 1, 0)))

#Construimos nuestra tabla de predicciones y datos
df <- cbind(df,
 data.frame(LinearSVM = ifelse(predict(linear.svm.fit) > 0, 1, 0),
 PolynomialSVM = ifelse(predict(polynomial.svm.fit) > 0, 1, 0),
 RadialSVM = ifelse(predict(radial.svm.fit) > 0, 1, 0),
 SigmoidSVM = ifelse(predict(sigmoid.svm.fit) > 0, 1, 0)))

#Reordenamos nuestros datos
predictions <- melt(df, id.vars = c('X', 'Y'))

#Graficamos los cálculos
ggplot(predictions, aes(x = X, y = Y, color = factor(value))) +
 geom_point() +
 facet_grid(variable ~ .)

 Ejemplo_2

Anteriormente mostré como al evaluar o clasificar correos el método SVM muestra un buen comportamiento al compararlo con el método clásico  Naive Bayes para analizar spam/no-spam.

Pero Naive Bayes no es el único método de clasificación entre los algoritmos de Machine Learning, así que hago otro comparativo entre diversas técnicas con datos de 3249 correos detectados como spam y se toman para evaluar los métodos 410 palabras contenidas en ellos. En el ejemplo anterior la base de los correos se encontraba ordenada y procesada, en este caso de hace uso de los textos para detectar y clasificar los correos.

Los datos de los correos se encuentran en formato de una matriz, de la cual se toma el 50% de datos para entrenamiento ( training set) y  prueba (test set) el otro 50%.

La cuestión ahora es elegir entre los otros métodos el que resulte mejor como modelo, para eso se hace uso de la validación cruzada y la regularización. Entonces se una el parámetro lambda y el MSE para determinar el lambda apropiado para el modelo logit.

#Cargamos los datos
load(file.path('data', 'dtm.RData'))

set.seed(1)

training.indices <- sort(sample(1:nrow(dtm), round(0.5 * nrow(dtm))))
test.indices <- which(! 1:nrow(dtm) %in% training.indices)

#Definimos los conjuntos de entrenamiento

train.x <- dtm[training.indices, 3:ncol(dtm)]
train.y <- dtm[training.indices, 1]

#Definimos los conjuntos de prueba

test.x <- dtm[test.indices, 3:ncol(dtm)]
test.y <- dtm[test.indices, 1]

#Borramos la matriz original
rm(dtm)

#Cargamos la librería glmnet para calcular logit
library('glmnet')

#Estimamos logit
regularized.logit.fit <- glmnet(train.x, train.y, family = c('binomial'))

lambdas <- regularized.logit.fit$lambda

performance <- data.frame()

#Extraemos los valores de lambda

for (lambda in lambdas)
{
 predictions <- predict(regularized.logit.fit, test.x, s = lambda)
 predictions <- as.numeric(predictions > 0)
 mse <- mean(predictions != test.y)
 performance <- rbind(performance, data.frame(Lambda = lambda, MSE = mse))
}

#Graficamos los resultados

ggplot(performance, aes(x = Lambda, y = MSE)) +
 geom_point() +
 scale_x_log10()

 La gráfica que se  obtiene es la siguiente:

Logit_Lambda

 

Ya que se elige el mejor modelo Logit se estima como se comportó al clasificar los correos.

#Resultados de Logit

best.lambda <- with(performance, max(Lambda[which(MSE == min(MSE))]))

mse <- with(subset(performance, Lambda == best.lambda), MSE)

mse
#[1] 0.06830769

Se observa que aproximadamente solo el 6% de los correos fueron mal clasificados.

Ahora calculo el modelo SVM para el kernel lineal y radial.

#Calculo de SVM
#Cargamos la librería que se usará
library('e1071')

#Estimamos le modelo para el kernel lineal
linear.svm.fit <- svm(train.x, train.y, kernel = 'linear')

predictions <- predict(linear.svm.fit, test.x)
predictions <- as.numeric(predictions > 0)

mse <- mean(predictions != test.y)

#Calculamos el valor del MSE
mse
#0.128

#Estimamos el modelo para el Kernel radial
radial.svm.fit <- svm(train.x, train.y, kernel = 'radial')

predictions <- predict(radial.svm.fit, test.x)
predictions <- as.numeric(predictions > 0)

#Calculamos MSE
mse <- mean(predictions != test.y)

mse
#[1] 0.1421538

Los dos modelos SVM tienen aproximadamente el 12% y 14% respectivamente de correos mal clasificados. Lo cual es mayor al error de clasificación con el modelo Logit.

Como último ejemplo se hace uso del algoritmo k-vecinos cercanos (k-nn), con k=50.

#Calculamos para 50 vecinos cercanos
#Usamos la librería class para la función knn
library('class')

#Estimamos el modelo
knn.fit <- knn(train.x, test.x, train.y, k = 50)

predictions <- as.numeric(as.character(knn.fit))

#Calculamos MSE
mse <- mean(predictions != test.y)

mse
#[1] 0.1396923

Se observa que este último modelo clasifica incorrectamente aproximadamente el 13 % de los correos.

Para mejorar el uso del modelo k-nn se hace la estimación del modelo para k desde 5 hasta 50 y  se elige el mejor modelo, considerando el modelo que tenga el  MSE con menor valor.

#Calculamos el k-nn modelo para varios k

performance <- data.frame()

for (k in seq(5, 50, by = 5))
{
 knn.fit <- knn(train.x, test.x, train.y, k = k)
 predictions <- as.numeric(as.character(knn.fit))
 mse <- mean(predictions != test.y)
 performance <- rbind(performance, data.frame(K = k, MSE = mse))
}

#Elegimos el mejor k
best.k <- with(performance, K[which(MSE == min(MSE))])

best.k
# 5
#Calculamos MSE para el mejor k
best.mse <- with(subset(performance, K == best.k), MSE)

best.mse
#[1] 0.09169231

Se obtiene que el mejor modelo de la familia k-nn resulta ser para k=5 y aproximadamente hace 9% de malas clasificaciones.

Usando como indicador de la calidad de los modelos al MSE se concluye que el modelo Logit tiene el 6% de errores. Así que el comportamiento de esta muestra de correos parece ser mejor modelado por técnicas lineales que no-lineales.

Lo que se observa con el ejemplo es que uno debe de realizar varios modelos para el mismo conjunto de datos y debe de tener algún criterio o indicador para determinar cual es el mejor modelo de entre todos los realizados. No debe uno de creer que una familia de modelos es mejor que la otra, solo alguna funciona mejor para ciertos datos que otra y nada más.

Siempre es conveniente tener más de un modelo y en caso de que dos puedan ser adecuados, es valido hacer ligeras variaciones de esos para determinar cual se adecua mejor al tipo de datos o fenómeno analizado.

Como ejemplo del tipo de gráficas que se obtiene de la librería kernlab hago un ejemplo. Cabe mencionar que esta librería y e1071 cuentan con varias funciones que concentrar algoritmos usuales en Machine Learning. Lo que las distingue es que la primera, kernlab; hace un concentrado los algoritmos dependientes de algún tipo de kernel.

Para el ejemplo se crean 300 puntos y se clasifican en dos grupos, análogo a los primeros ejemplos se busca la curva que separa en dos clases a los datos.

#Generamos los datos

n<-300 #Número de puntos
p<-2 #Dimensiones

sigma<-1 #Varianza de la distribución
meanpos<-0 #Centro de la distribución de los puntos positivos
meanneg<-3 #Centro de la distribución de los putnos negativos

npos<-round(n/2) #Número del ejemplo positivos
nneg<-n-npos #Números negativos del ejemplo

#Generación de los postivos y negativos
xpos<-matrix(rnorm(npos*p, mean=meanpos,sd=sigma),npos,p)
xneg<-matrix(rnorm(nneg*p, mean=meanneg,sd=sigma),nneg,p)
x<-rbind(xpos,xneg)

#Generación de etiquetas

y<-matrix(c(rep(1,npos),rep(-1,nneg)))

#Visualización de datos
plot(x,col=ifelse(y>0,1,2))
legend("topleft",c('Positivos','Negativos'),col=seq(2),pch=1,text.col=seq(2))

Ejemplo_3 La gráfica muestra la mezcla de puntos en color rojo y negro. Primero se estima el modelo lineal de SVM para determinar la clasificación,  en este caso se toma un subconjunto de entrenamiento y otro de prueba.

#Separamos los datos en entrenamiento y prueba

ntrain <- round(n*0.8) #Entrenamiento
tindex <- sample(n,ntrain) #Muestra de índices para entrenamiento
xtrain <- x[tindex,]
xtest <- x[-tindex,]
ytrain <- y[tindex]
ytest <- y[-tindex]
istrain=rep(0,n)
istrain[tindex]=1 

#Visualización

plot(x,col=ifelse(y>0,1,2),pch=ifelse(istrain==1,1,2))
legend("topleft",c('Positive Train','Positive Test','Negative Train','Negative Test'),col=c(1,1,2,2),pch=c(1,2,1,2),text.col=c(1,1,2,2))

library(kernlab)

#Entrenamiento de SVM

svp <- ksvm(xtrain,ytrain,type="C-svc",kernel='vanilladot',C=100,scaled=c())

#Graficación

plot(svp,data=xtrain)

 Ejemplo_3_1

Ejemplo_3_2

 

En la primera gráfica se muestra como se comportan los subconjunto de entrenamiento y prueba. En la segunda gráfica se muestra como el método SVM separa los datos y se queda con los negritos dentro de lo que se conoce como margen. Las clases en las cuales separa los datos les da color rojo y azul. En general cada librería cuenta con algunos datos los cuales se encuentran procesados y son un recurso para probar con las funciones y métodos.

Un ejemplo para usar SVM para regresión

Para este ejemplo uso datos de R Project correspondientes a un indicador financiero Alemán; DAX, estos corresponden a los cierres entre 1991 y 1998.

#Se cargan los datos
data("EuStockMarkets")
#Exploramos los datos
head(EuStockMarkets)

#Se toma solo solo el indicador DAX y se construye 
#una variable para el tiempo, t

t=1:length(EuStockMarkets[,1])
M=data.frame(EuStockMarkets[,1],t)
colnames(M)<-c("DAX","Tiempo")

#Se construye el modelo con SVM
Modelo=ksvm(M$Tiempo,M$DAX)

#Se hace la gráfica del modelo
plot(M$Tiempo,M$DAX,type="l", main="Ejemplo de SVM para regresión", xlab="Tiempo",ylab="DAX")
lines(M$Tiempo,predict(Modelo,M$Tiempo),col="red")


SVM_Regresion_BlogLa imagen muestra el resultado que predice un modelo regresivo con la familia de técnicas de SVM, este modelo no es el mejor, ya que se tendría que explorar el comportamiento del indicador DAX y hacer algunas pruebas para elegir el mejor kernel para este modelo. Pero la intención es mostrar que SVM también puede ser usado para construir modelos regresivos, la desventaja ante modelos de regresión es la expresión o ecuación, pero como mencioné anteriormente esto realmente es relevante dependiendo de los datos  y el origen de la necesidad de modelar estos datos.

Espero sirva de guía lo mostrado en la entrada, en las referencias se encuentran muchos detalles que no menciono.

Referencias:

1.-http://shop.oreilly.com/product/0636920018483.do

2.-http://www.amazon.com/Machine-Learning-Probabilistic-Perspective-Computation/dp/0262018020

3.-http://research.microsoft.com/en-us/um/people/cburges/papers/svmtutorial.pdf 

4.-http://www.springer.com/gp/book/9780387987804

5.-http://www.kernel-machines.org/

6.-http://www.support-vector-machines.org/

7.-http://www.csie.ntu.edu.tw/~cjlin/libsvm/

8.-El capítulo 7, http://www.springer.com/gp/book/9780387310732

9.-El capitulo 12, http://statweb.stanford.edu/~tibs/ElemStatLearn/

Optimización

Sobre Optimización

Para entender en general las ideas de esta entrada basta con tener algunos conocimientos de Calculo Diferencial, de una y varias variables. Si no se cuenta con ideas o conocimientos de ello, no implica el no entender los algoritmos de optimización que se comentan o por lo menos visualizar como funcionan.

Parto de ideas sencillas hasta dar unos ejemplos de optimización estocástica. Espero poder trasmitir la idea de los métodos y que puedan servirles de ayuda para investigar más respecto al tema.

Omito hablar sobre el proceso ” gradiente descendiente estocástico” pero dejo referencias donde pueden consultar en qué consiste este método y por que es que se hace amplio uso en técnicas de Machine Learning.

 La idea de optimizar  seguro es familiar cuando vemos como nuestras madres resuelven la siguiente pregunta: ¿cómo hacer para comer carne con solo cierta cantidad de dinero?

Aunque parece broma la pregunta anterior es una optimización, tenemos que encontrar las cantidades de bienes que podemos comprar considerando que solo tenemos una cantidad de dinero y debemos de incluir en la lista de bienes la carne. Obviamente esto no dice que nuestras madres hagan uso de las técnicas que comento esta entrada, pero espero que de una idea burda de lo que es optimizar.

Otro ejemplo es cómo maximizar las ganancias de una empresa o negocio, considerando R cantidad de ingresos y C cantidad de costos. Claramente este problema implica conocer qué genera nuestros ingresos y costos, para plantear un modelo idealizado sobre la empresa.

En general la idea de optimizar es encontrar los máximos o mínimos de una función. Para tener una idea de qué significa hago dos gráficas como ejemplo en R project.

#Creamos una función
#Ejemplo 1-Mínimo de una función

Func_1<-function(x){
 x[1]^2 
      }

resultado<-optimize(Func_1,c(-5,5),maximum=FALSE,tol=1e-8) 

#Vemos los resultados
print(resultado$minimum)
print(resultado$objective)

#Generamos el rango de valores para x y y=f(x)
x<-seq(-5,5,length.out=100)
y<-Func_1(expand.grid(x))

#Graficamos identificando el mínimo

plot(x,y,xlab="x",ylab="f(x)",col="4")
points(resultado$minimum,resultado$objective, col="red",pch=19)
rect(resultado$minimum-0.3, resultado$objective-0.7, resultado$minimum+0.3,resultado$objective+0.7)

Ejemplo_opt_mínimo

La gráfica muestra en color rojo el mínimo de la función, lo obtenemos por medio de la función optimize. La cual necesita una función objetivo, un rango de valores para optimizar y un nivel de tolerancia en el cálculo. Esto último es el número que utilizará para realizar su aproximación. Recomiendo revisar la información de la función por medio del comando help.

Cuando pedimos el resultado por medio de print(resultado$minimun)en el código, se observa que el valor es 0, lo cual es evidente al ver la gráfica. Hago otro ejemplo con una función diferente.

#Ejemplo 2
#Creamos una función
Func_2<-function(x){
 x^2+6*sin(x*pi*2)
 }
 
#Optimizamos
resultado2_1<-optimize(Func_2,c(-5,5),lower="-1",upper="3",maximum=FALSE,tol=.Machine$double.eps^0.25) 
resultado2<-optimize(Func_2,c(-5,5),maximum=FALSE,tol=1e-8)
#Vemos los resultados
print(resultado2$minimum)
print(resultado2$objetive)

#Rango de calores para nuestras gráficas 
x2<-seq(-5,5,length.out=100)
y2<-Func_2(expand.grid(x2)) 

#Graficamos
plot(x2,as.matrix(y2)[,1],xlab="x",ylab="f(x)",col="4")
points(resultado2$minimum,resultado2$objective, col="red",pch=19)
rect(resultado2$minimum-0.3, resultado2$objective-0.7, resultado2$minimum+0.3,resultado2$objective+0.7)

points(resultado2_1$minimum,resultado2_1$objective,col="green",pch=19)
rect(resultado2_1$minimum-0.3,resultado2_1$objective-0.7,resultado2_1$minimum+0.3,resultado2_1$objective+0.7)

Ejemplo_opt_mínimo_2

La diferencia de esta última gráfica es que se tienen dos puntos, rojo y verde, que a primera vista parece que el verde es el “mínimo” de toda la gráfica. Pero lo que realmente hice fue acotar el cálculo del optimización, por lo cual si observamos en la función  optimize hice distintas cotas (lower, upper) para que se ejecutara la función. Entonces lo que se ve son mínimos locales y dependen de la restricción que se haga a los cálculos.

¿Qué es lo que hace la función optimize? Da el punto máximo o mínimo según las restricción para el cálculo o en otro caso regresa el máximo o mínimo global.

La diferencia en local y global, es que el primero se restringe a una pequeña región, ejemplo a las x que están entre -2 y -1, el punto rojo de la gráfica es un mínimo local. Por otro lado, global se refiere a que para todo el rango de las x’s que se revisan no se puede encontrar otro valor de f(x) menor. En el ejemplo el punto verde de la segunda gráfica es un mínimo global y el punto rojo de la primera también.

Los dos ejemplos muestran gráficamente cual es la idea de obtener un mínimo que depende de una sola variable x. Pero en la mayoría de problemas de optimización no son funciones de una sola variable, sino de varias variables. Ejemplo, cuando  se busca maximizar ganancias o minimizar costos, es difícil que dependa de un solo producto que se vende o el costo de un solo insumo.

Entonces, la idea de optimizar se puede generalizar para varias variables. Pero en este caso se tienen varios métodos para hacer el cálculo y cada uno tienen ventajas y desventajas, depende del tipo de problema o función que  se desea optimizar.

¿Qué tiene que ver todo esto con Machine Learning?

 En general todos los algoritmos de Machine Learning tienen tres piezas, datos de entrada, datos de salida y una función entre los datos. Pero por medio de un conjunto de datos de entrenamiento y un algoritmo de aprendizaje se trata de obtener la mejor función entre los datos.

Entonces, la optimización está ligada a los algoritmos, ya que siempre en los algoritmos  se tiene una función objetivo y un conjunto de parámetros que se busca estimar. Por lo cual es fundamental la optimización al aplicar los algoritmos de Machine Learning.

Por ejemplo, en la regresión lineal. Sobre la cual en otra entrada se mencionó que la función lm permite estimarla. Se toman datos de altura y peso de personas, los datos se pueden descargar desde de Github.

Al usar la función lm  no se muestra lo que hace o como estima los valores de la regresión  a y b, donde la ecuación es:

Peso=a+b*Altura

Lo que se puede hacer es estimar los valores de los parámetros por medio de una optimización y comparar con los valores que se estiman en R project por medio de la función lm.

#Ejemplo de estimación de Parámetros
#Definimos la función de la regresión

Altura.a.peso <- function(altura, a, b)
{
 return(a + b * altura)
}

#Función de error cuadrado
error.cuadrado <- function(altura.peso, a, b)
{
 predic <- with(altura.peso, Altura.a.peso(Altura, a, b))
 errors <- with(altura.peso, Peso - predic)
 return(sum(errors ^ 2))
}

#Optimización

optim(c(0, 0),function (x){error.cuadrado(altura.peso, x[1], x[2])})

El código anterior hace dos cosas. La primera es definir la función de la cual se busca estimar sus parámetros a y b, pero esa función solo define rectas y para elegir la que mejor describe la relación entre  los datos se necesita otra función, error.cuadrado.

La función error.cuadrado como se observa depende de la función de Altura.a.peso. Así que se tiene solo una función compuesta con otra y  se busca en qué valores se minimiza esa función.

La función error.cuadrado permite encontrar la mejor recta entre los valores de Peso contra Altura. Además indica que se encontró la recta que tienen el menor error cuadrado medio, que es medido entre el Peso menos el valor que estimamos mediante la función Altura.a.peso.

Comparo los valores que se encuentran con la estimación por optimización y por la función lm.

#Comparación entre las estimaciones

 >lm(Peso~Altura,data=altura.peso)

Call:
lm(formula = Peso ~ Altura, data = altura.peso)

Coefficients:
(Intercept) Altura 
 -350.737 7.717 


> optim(c(0, 0),function (x){error.cuadrado(altura.peso, x[1], x[2])})
$par
[1] -350.786736 7.718158

$value
[1] 1492936

$counts
function gradient 
 111 NA 

$convergence
[1] 0

$message
NULL

Se nota que los valores de la regresión por medio de lm y los valores estimados por la función optim son prácticamente los mismos. Entonces se obtiene mediante optimización de una función objetivo error.cuadrado  el valor de sus parámetros a, b.

Se observa que la optimización fue hecha para una función de dos variables, el método que se usa para optimizar es Nelder-Mead, más adelante hablo de él con algo de detalle.

En conclusión, lo que parecía que operaba como una caja-negra cuando se usa la función lm, ahora se puede entender o tener una visión de cómo posiblemente es el proceso de estimación de los parámetros.

Aprovechando el ejemplo de la regresión lineal se puede modificar la función de error para estimar lo que se llama Rigde Regression. Es otro tipo de regresión en la cual para estimar los parámetros se afecta su error cuadrado o se penalizan los errores. El objetivo de esta técnica es evitar la sobre estimación de los datos y disminuir la varianza del ajuste lineal.

La Ridge Regression implica agregar un nuevo parámetro, lambda; el cual para hacer una buena elección se requiere hacer uso de la Regularización y Cross-Validation[1]. Se haciendo la modificación a la estimación para revisar el comportamiento al optimizar.

Haciendo dos estimaciones de Ridge Regression y se compara la estimación con respecto a las que se realizo anteriormente.

#Calculo de Ridge Regression

ridge.error <- function(altura.peso, a, b, lambda)
{
 predic <- with(altura.peso, Altura.a.peso(Altura, a, b))
 errores <- with(altura.peso, Peso - predic)
 return(sum(errores ^ 2) + lambda * (a ^ 2 + b ^ 2))
}

lambda=1
optim(c(0, 0),function (x){ridge.error(altura.peso, x[1], x[2], lambda)})

$par
[1] -340.434108 7.562524
lambda=0.75
optim(c(0, 0),function (x){ridge.error(altura.peso, x[1], x[2], lambda)})

$par
[1] -343.13351 7.60324

El ejemplo muestra el uso de la optimización al estimar los parámetros de un modelo. En la Rigde Regression se definió el valor de lambda, lo cual hace que se tenga ese parámetro “libre”  en el modelo y por ello es necesario para un buen uso de esta técnicas elegir la mejor lambda por medio de regularización y cross-validation.

Entonces al usar el paquete ridge de R project para estimar la ridge regression, ya no parece del todo extraño como hace el cálculo el software. Podría usar la función linearRidge del paquete ridge y comparar por medio de optimización, regularización y cross-validacón la elección del lambda.

Algo que no hice pero que puede resultar ser un buen ejercicio es revisar el artículo de investigación de Erika Cule y estudiar el método de elección de lambda contra otro método para elegirla.

Ejemplo_opt_regresion

La gráfica muestra como cambia ligeramente la pendiente de las rectas. Haciendo solamente una modificación en la estimación de los errores, regresión lineal y ridge regression.

Espero que este ejemplo haga notar que los algoritmos de Machine Learning  que parecen cajas negras, en el fondo están combinando cosas, y una de esas es un proceso de optimización para estimar los parámetros.

Este tema un trabajo serio al momento de procesar volúmenes masivos de datos, quienes tengan interés pueden revisar las investigaciones del grupo de optimización de la universidad Carnegie-Mellon.

Cabe mencionar que en las investigaciones sobre optimización su relevancia puede implicar una disminución en el costo computacional o en la mejora de algoritmos.

Ejemplos de técnicas de Optimización

Lo siguiente que comparto en esta entrada son técnicas de optimización convexa y dos técnicas de optimización estocástica. Para algunas solo pondré ejemplos de cómo se estima y como interpretar el resultado siguiendo unos ejemplos sencillos publicados por Jason Brownlee en su página. En otros pondré ejemplos de como usarlas de manera más detallada.

La explicación de cada uno de los métodos no es exhaustiva, tanto por que está fuera de mi alcance como de la intensión de esta entrada.

Dejo referencias en los comentarios y al final de la entrada algunas ligas a notas, reportes de investigación y libros para consultar más sobre el tema.

Método 1.-Búsqueda de la sección dorada

No repito el código, pero la primera gráfica que se hizo en esta entrada hacía uso de la función optimize y muestra un punto rojo que es el mínimo de la función. El método por medio del cual estima ese punto es conocido como  búsqueda de la sección dorada.

Su uso se limita a función que dependen de una sola variable y se busca el punto óptimo entre dos puntos que se definen del rango de las x’s. El nombre provienen de la relación entre los tres primeros puntos que toma el método, ya que esto tienen una proporción con el número conocido como razón dorada o número áureo.

Entre los limitantes de dicha técnica es que solo funciona para funciones de una sola variable, de preferencia que sea una función convexa y que se pueda identificar que se tienen un único óptimo entre la región que define para realizar el cálculo.

Función Rosenbrock

Para tres de los siguientes ejemplos se hace uso de la función Rosenbrock. Es una función no convexa que se usa para probar los métodos de optimización o rendimiento de los algoritmos. En general se puede usar de dos variables o se puede generalizar para n-variables. Forma parte de un conjunto de funciones que son usadas para hacer pruebas con técnicas de optimización.

La ecuación en general de la función para dos variables es:

f(x,y)=(a-x)^2+b(y-x^2)^2

La visualización de la función para un par de valores a=1 y b=100.

#Función Rosenbrock
#La función es f(x,y)=(1-x)^2+100*(y-x^2)^2

nx<-21
ny<-21

x<-seq(-5,5,length=nx)
y<-seq(-5,5,length=ny)

#Se genera la salida de los valores en z

z<-outer(x,y,function(x,y)(1-x)^2+100*(y-x^2)^2)

#Escala de colores
hgt <- 0.25 * (z[-nx,-ny] + z[-1,-ny] + z[-nx,-1] + z[-1,-1])
hgt <- (hgt - min(hgt))/ (max(hgt) - min(hgt))

#Se construye la gráfica en perspectiva con cierto ángulo para visualizarla
persp(x, y, z, col =cm.colors(10)[floor(9* hgt + 1)],theta = 35,phi=30)

Función_Rosenbrock

Método 2.-Nelder-Mead

El método es adecuado para funciones no lineal de varias variables y se busca el óptimo sin restricciones en el dominio. La herramienta que usa para buscar el óptimo en su algoritmo es una estructura geométrica conocida como simplex. El método para una función de dos variables parte de la elección de 3 puntos iniciales de los cuales bajo un proceso de varias etapas se va modificando el simplex hasta aproximarse al punto óptimo de la función. El proceso de transformación del simplex es la parte principal del algoritmo.

Es adecuado para funciones de varias variables y al no requerir la derivada para su algoritmo puede ser usada para funciones no diferenciables, no continuas o con ruido. Un resumen de sus ventajas y desventajas se pueden consultar en, Nelder-Mead_algorithm.

#Creamos una función

#Definición de una función Rosenbrock, el optimo de ubica en (1,1)

rosenbrock <- function(v){ 
 #Valores de a=1, b=100
(1 - v[1])^2 + 100 * (v[2] - v[1]*v[1])^2
}

# Identificación del mínimo local por el método Nelder-Mead 
resultado3 <- optim(c(runif(1,-3,3), runif(1,-3,3)),rosenbrock,NULL,method="Nelder-Mead",control=c( maxit=100,reltol=1e-8,alpha=1.0, beta=0.5,gamma=2.0))

#Revisión de resultados
print(resultado3$par) 
print(resultado3$value) 
print(resultado3$counts)

#Gráfica de cortorno
x3<- seq(-3, 3, length.out=100)
y3<- seq(-3, 3, length.out=100)
z3<- rosenbrock(expand.grid(x3, y3))
contour(x3, y3, matrix(log10(z3), length(x3)), xlab="Valores de X",ylab="Valores de Y")
#Punto óptimo
points(resultado3$par[1], resultado3$par[2], col="red", pch=19)
#Úbicación del punto óptimo
rect(resultado3$par[1]-0.2, resultado3$par[2]-0.2, resultado3$par[1]+0.2, resultado3$par[2]+0.2, lwd=2)

Ejemplo_opt_Nelder-Mead2

Método 3.-Gradiente Descendente o Algoritmo de Descenso Rápido.

El métodos hace uso del gradiente de una función y cuando busca el mínimo hace uso del negativo del gradiente para aproximarse al óptimo.

Funciona bien si la función es convexa, ya que el método estimará el mínimo global. El método es muy sensible a la elección del punto inicial para hacer la estimación y como depende de un parámetro para ir haciendo los avances, es susceptible a malas elecciones de él.

Como referencias dejo las siguientes ligas: Aplicación del método a procesamiento de imágenes y curso de Andrew Ng en línea.

Usamos el método para calcular el mínimo de una función que tienen su mínimo global en el punto (0,0), el cual se visualiza en la imagen.

Ejemplo_función_min_global

#Creamos una función
#Se define una función, el óptimo se ubica en (0,0)

Func_3<- function(x) {
x[1]^2 + x[2]^2
}

#Definición de la derivadas
derivada <- function(x) {
c(2*x[1], 2*x[2])
}

#Definición del método del gradiente
gradient_descent <- function(func, derv, start, step=0.05, tol=1e-8) {
 pt1 <- start
 grdnt <- derv(pt1)
 pt2 <- c(pt1[1] - step*grdnt[1], pt1[2] - step*grdnt[2])
 
 while (abs(func(pt1)-func(pt2)) > tol) {
 pt1 <- pt2
 grdnt <- derv(pt1)
 pt2 <- c(pt1[1] - step*grdnt[1], pt1[2] - step*grdnt[2])
 print(func(pt2)) #regresa el valor en el punto
 }
 pt2 # Regresa el valor óptimizado
 }

#Método del Gradiente Descendente
resultados4 <- gradient_descent(Func_3, derivada,c(runif(1,-3,3), runif(1,-3,3)), 0.05,1e-8)

#Resumen de resultados
print(resultados4) 
print(Func_3(resultados4)) 

#Gráfica de contorno
x4 <- seq(-3, 3, length.out=100)
y4 <- seq(-3, 3, length.out=100)
z4 <- Func_3(expand.grid(x4, y4))
contour(x4, y4, matrix(z4, length(x4)), xlab="x",ylab="y")
#Dibujamos el punto óptimo
points(resultados4[1], resultados4[2], col="red", pch=19)
#Dibujamos el cuadrado del óptimo
rect(resultados4[1]-0.2, resultados4[2]-0.2, resultados4[1]+0.2, resultados4[2]+0.2, lwd=2)

Ejemplo_gradiente_descendente

Método 4.-Gradiente Conjugado

Es adecuado para funciones no lineales de varias variables. Le idea del método es construir una base ortogonal de vectores y hacer uso de ella para hacer la búsqueda del óptimo de un modo más eficiente. Aspectos teóricos de manera breve puede consultar en las notas en español del Prof. Maximenko y una buena referencia para ver más detalles recomiendo las notas del Prof.Shewchuk.

Igual que otros métodos su eficiencia se ve mermada ante una mala elección del punto inicial en  funciones no  convexas, pero una de sus ventajas es que no hace uso de la matriz Hessiana en su algoritmo lo cual permite hacer uso del método para funciones de muchas variables de manera eficiente.

Para hacer el ejemplo hacemos uso de la función Rosenbrock.

#Definimos la función Rosenbrock

rosenbrock <- function(v) { 
(1 - v[1])^2 + 100 * (v[2] - v[1]*v[1])^2
}

#Definición de la derivada de la función
derivada2 <- function(v) {
c(-400 * v[1] * (v[2] - v[1]*v[1]) - 2 * (1 - v[1]), 
200 * (v[2] - v[1]*v[1]))
}


Ejemplo_opt_gradente_conjugado2

Método 5.-BFGS

Este método hace uso tanto del gradiente como de una aproximación a la inversa de la matriz Hessiana de la función, esto para hacer una aproximación al cálculo de la segunda derivada. Por ser un método de aproximación a la segunda derivada se dice que es un método Quasi-Newtoniano. El problema con el método es el costo computacional para funciones de muchas variables, ya que se almacena una matriz cuadrada de datos tan grande como el cuadrado de la cantidad de variables.

Es adecuado para funciones no lineales de varias variables y la búsqueda del óptimo sin restricciones. Si bien el método hace que la convergencia al óptimo sea rápida por ser de aproximación a la segunda derivada, el costo de procesamiento es bastante alto.

Se pueden consultas las notas de Galen Andrew o de Peter Blomger, para ver detalles del método y el algoritmo.

#Método BFGS

#Definición la función Rosenbrock, se óptimiza en (1,1)
rosenbrock <- function(v) { 
(1 - v[1])^2 + 100 * (v[2] - v[1]*v[1])^2
}
 
#Definición del Gradiente
derivada4 <- function(v) {
c(-400 * v[1] * (v[2] - v[1]*v[1]) - 2 * (1 - v[1]), 200 * (v[2] - v[1]*v[1]))
}
 
#Calculo del método BFGS 
resultados6 <- optim(c(runif(1,-3,3), runif(1,-3,3)),rosenbrock, derivada4,method="BFGS",control=c(maxit=100,reltol=1e-8))
 
#Resumen de resultados
print(resultados6$par)
print(resultados6$value) 
print(resultados6$counts)
 
#Gráfica de contorno
x5 <- seq(-3, 3, length.out=100)
y5 <- seq(-3, 3, length.out=100)
z5 <- rosenbrock(expand.grid(x5, y5))
contour(x5, y5, matrix(log10(z5), length(x5)), xlab="x", ylab="y")

points(resultados6$par[1], resultados6$par[2], col="red", pch=19)
# draw a square around the optima to highlight it
rect(resultados6$par[1]-0.2, resultados6$par[2]-0.2, resultados6$par[1]+0.2, resultados6$par[2]+0.2, lwd=2)

Ejemplo_BFGS

Los métodos anteriores son considerados de optimización convexa, para ver detalles teóricos recomiendo consultar un texto especializado en el tema. Una buena referencia es el libro de Stephen P. Boyd y Lieven Vandenberghe el cual se puede descargar libremente desde su página. El segundo autor tienen dos cursos en línea con formato de presentaciones; optimización convexa y métodos de optimización para sistemas de larga escala.

Optimización Estocástica (o ¿computación evolutiva?)

Los siguientes ejemplos que comparto son abordados con técnicas de optimización estocástica, las cuales también son estudiadas en computación evolutiva. Por supuesto que no doy explicaciones exhaustivas, pero dejo referencias al final de la entrada.

Los métodos estocásticos que generalmente se usan son: algoritmos genéticos (genetic algorithm)  y recorridos simulados (Simulated anneling).

En el ejemplo considero el problema siguiente, si se  tiene un grupo de viajeros y todos parten de lugares distintos pero se dirigen al mismo sitio, si todos están de acuerdo en reunirse en un punto y se esperaran para compartir transporte terrestre, esto implica que entre los posibles vuelos tienen que elegir aquellos que mejor les convenga en costo, tiempo de espera, horarios de llegada y salida. Así que se tienen varios factores a considerar en la optimización.

Como cualquier problema de optimización se requiere tener la función objetivo, en este ejemplo se le llamará función de costo.

Los datos pueden ser descargados desde la página de Toby Segaran, cuenta con con una versión del código en python. Hice una adaptación en R y considerando las restricciones del problema y de los datos.

Los pasos para abordar el problema son los siguientes:

  1. Se representarán las soluciones como cadenas de números.
  2. Se presentaran las soluciones como un tabla.
  3. Se definirán una función objetivo o función de costo.
  4. Se estimaran soluciones por medio del método búsqueda aleatoria, recorrido simulados y algoritmos genéticos.

Paso 1

Se representa las solución como una cadena de 12 dígitos que pueden estar entre 1 y 9. Ejemplo, solución=[1,3,4,5,8,10,2,3,7,9,3,5], representa que la persona uno toma el primero vuelo y regresa en su tercer vuelo posible. La persona dos toma el cuarto vuelo y regresa en su quinto vuelo posible, así sucesivamente para las 6 personas.

#Carga de los datos desde el directorio
viajes.path<-file.path("data","schedule.txt")

viajes<-read.table(viajes.path,sep=",",stringsAsFactor=FALSE)

colnames(viajes)<-c('Origen','Destino','Hora de Salida','Hora de arrivo','Precio')

#Generamos lista de viajeros

gente=list(c('Seymour','Franny','Zooey','Walt','Buddy','Les'),c('BOS','DAL','CAK','MIA','ORD','OMA'))

#Destino al que viajan
Destino=c('LGA')

#Definimos una función para medir el tiempo, será necesaria en otras funciones
#Librería para manipulas fechas
library(lubridate)

#Función que nos regresa los minutos
minutos<-function(horario){
 min<-hour(hm(horario))*60+minute(hm(horario))
 return (min)}

Paso 2

Teniendo los datos cargados, se procede a representar por una pequeña tabla la solución al problema. Para ello se define la función impri.horario.

#Impresión de soluciones

impri.horarios<-function(r){

#r es la cadena con datos de los viajeros y lugars
 for( i in 1:(length(r)/2)){
 nombre=gente[[1]][i]
 origen=gente[[2]][i]
 Destino=c('LGA')
 
 salida=viajes[which(viajes$Origen==origen),][r[i*2-1],]
 regreso=viajes[which(viajes$Origen==Destino),][r[i*2],]
 
 a<-c(nombre,origen,paste(salida[,3],"-",salida[,4],sep=""," ",salida[,5]),paste(regreso[,3],"-",regreso[,4],sep=""," ",regreso[,5]))
 print(a,sep="\t")
 
 } 
 
} 

Paso 3

Se define una función objetivo o la función de costo, es un ejemplo solamente, por lo cual la función puede ser modificada para una situación real de optimización. Pero cabe mencionar que la definición de ella es la parte más importante para aplicar estos algoritmos.

#Creamos una función a optimizar
#Función Objetivo o Función de Costo

costo.tiempo<-function(r){
 precio.total=0
 ultimo.arrivo=0
 salida.temprana=24*60
 
 for(i in 1:(length(r)/2)){
 #nombre=gente[[1]][i]
 origen=gente[[2]][i]
 Destino=c('LGA')
 salida=viajes[which(viajes$Origen==origen),][r[i*2-1],]
 regreso=viajes[which(viajes$Origen==Destino),][r[i*2],]
 precio.total=precio.total+salida[,5]+regreso[,5]
 
 if(ultimo.arrivo<minutos(salida[,4])){ultimo.arrivo=minutos(salida[,4])}
 if(salida.temprana>minutos(regreso[,3])){salida.temprana=minutos(regreso[,3])}
 
 }
 tiempo_espera=0 
 for( i in 1:(length(r)/2)){
 origen=gente[[2]][i]
 Destino=c('LGA')
 salida=viajes[which(viajes$Origen==origen),][r[i*2-1],]
 regreso=viajes[which(viajes$Origen==Destino),][r[i*2],]
 
 tiempo_espera=tiempo_espera+(ultimo.arrivo-minutos(salida[,4]))
 tiempo_espera=tiempo_espera+(minutos(regreso[,3])-salida.temprana)
 }
 if(ultimo.arrivo>salida.temprana){precio.total=precio.total+50} 
 
 return(precio.total+tiempo_espera)
 }

Si se hace un prueba con esta función se nota que devuelve un número entero, el cual representa el valor de la solución o configuración que asignamos a los viajantes.

#Ejemplo del uso de la función de costo
#Definimos una posible solución
 vector
 [1] 6 7 4 2 6 7 4 2 6 2 6 8
> costo.tiempo(vector)
[1] 3806

Paso 4

Para comparar métodos uso el método de búsqueda aleatoria. Consiste básicamente en tomar posibles soluciones y comparar entre las que se van encontrado, así con un proceso de comparación se elige la solución que minimiza  la función objetivo. Como es de esperar no es un método eficiente, pero permite acercarse a las ideas de la optimización estocástica.

#Búsqueda Aleatoria

optimizacion.aleatoria<-function(funcost){
 #dominio va ser igual al número máximo que puede tomar cada entrada
 best=999999999
 bestr=list()
 for( i in (0:1000)){
 r=sample(10,size=12,replace=TRUE) 
 cost=funcost(r)
 
 if(cost<best){
 best=cost
 bestr=r}
 }

 return(bestr)}

Cabe mencionar que el método fue adaptada para el problema, el algoritmo es estándar. Se ponemos a prueba el método y se evalúa la función de costo.

#Resultados del método Búsqueda Aleatoria

> optimizacion.aleatoria(costo.tiempo)
 [1] 8 9 6 7 2 8 8 7 6 8 7 8

> solucionba=c(8, 9, 6, 7, 2, 8, 8, 7, 6, 8, 7, 8)

> costo.tiempo(solucionba)
[1] 3689

> impri.horarios(solucionba)
[1] "Seymour" "BOS" "17:11-18:30 108" "18:25-20:34 205"
[1] "Franny" "DAL" "13:54-18:02 294" "15:07-17:21 129"
[1] "Zooey" "CAK" "8:27-10:45 139" "16:35-18:56 144"
[1] "Walt" "MIA" "17:07-20:04 291" "15:07-17:21 129"
[1] "Buddy" "ORD" "14:22-16:32 126" "16:35-18:56 144"
[1] "Les" "OMA" "15:03-16:42 135" "16:35-18:56 144"

Ahora usamos el algoritmo de recorrido simulados, tienen su origen en un problema de física. Por lo cual algunos de sus parámetros prevalecen en uso por analogía del fenómeno físico.

La idea es iniciar con una solución aleatoria y con un parámetro de temperatura alto, este último representa en la búsqueda la guía de un proceso de “enfriamiento” donde conforme decrece se van comparando los valores de la función objetivo. De los cuales se elige una posible solución al óptimo global que se espera obtener al tener un nivel de temperatura bajo.

El método es una adaptación al método Metropolis-Hastings, por lo cual en el algoritmo aparece una exponencial para comparar un valorar aleatorio y continuar en las comparaciones de los valores de la función objetivo.

#Algoritmo recorrido simulado
annealingoptimize<-function(costf,T=10000.0,cool=0.95,step=1){
 #Creamos una solución aleatória
 vec=sample(9,size=12,replace=TRUE)

 while(T>0.1){
 #Elejimos un índice
 i=sample(8,size=1)
 #Definimos la dirección del algortimos
 x=c(-1,0,1)
 dir=sample(x,1) 
 vecb=vec
 vecb[i]=vecb[i]+dir
 
 if(vecb[i]<0){vecb[i]=0}
 if(vecb[i]>9){vecb[i]=9}
 
 ea=costf(vec)
 eb=costf(vecb)
 
 
 if(eb<=ea){
 vec=vecb
 }
 else
 { 
 a=runif(1,min=0,max=1)
 
 if(a<exp((-eb-ea)/T))
 {
 vec=vecb
 
 } 
 else
 {
 vec=vec}
 
 
 }
 
 T=T*cool
 
 }
 return(vec)
 
 }

Hacemos una prueba de la optimización por este método. Los resultados fueron los siguientes.

#Método de Recorrido Simulado

annealingoptimize(costo.tiempo,T=10000.0,cool=0.95,step=1)
 [1] 8 7 4 7 8 3 7 3 1 5 8 7

> solucionrs=c(8,7,4,7,8,3,7,3,1,5,8,7)

#Evaluamos nuestra función de costo
> costo.tiempo(soucionrs)
[1] 4573
#Imprimimos nuestros vuelos

> impri.horarios(soucionrs)
[1] "Seymour" "BOS" "17:11-18:30 108" "15:07-17:21 129"
[1] "Franny" "DAL" "10:30-14:57 290" "15:07-17:21 129"
[1] "Zooey" "CAK" "17:08-19:08 262" "9:31-11:43 210" 
[1] "Walt" "MIA" "15:34-18:11 326" "9:31-11:43 210" 
[1] "Buddy" "ORD" "6:05-8:32 174" "12:31-14:02 234"
[1] "Les" "OMA" "16:51-19:09 147" "15:07-17:21 129"

Por último hago un pequeño código de un algoritmo genético, este fue adaptado para el problema que se está resolviendo. Lo que se debe de mencionar son las etapas que conforman dicho algoritmo: Mutación, Cruza, Generación de Población, Elitismo y procesamiento.

Cada etapa tienen varias opciones para ser realizada, la mutación puede ser solo una variación o afectación a la cadena de entrada. La cruza puede tener varios métodos para implementarse, y la idea base es “pegar” dos cadenas para obtener una nueva. La Población, tal cual como su nombre es, genera una cantidad de cadenas obtenidas de manera aleatoria y por “elitismo” elegir las mejores, es decir por un criterio se discrimina entre las cadenas.

Después de tener una Población inicial y haber elegido entre ellas,  se usa la mutación y la cruza para ir “mejorando” la población y por ultimo se  aproxima al óptimo.

El código que muestro es una versión simple y hasta cierto punto fea o mala, ya que hay pasos que pueden hacerse de manera más simplificada o compacta. Pero preferí mostrar todos los pasos por simples que fueran, pensando dar visibilidad en lugar de elegancia.

#Algoritmo Genético

genetic.optimize<-function(costf,popsize=50,step=1,mutprod=0.2,elite=0.2,maxiter=100){
 
 #Mutación
 mutacion<-function(vec){
 i=sample(12,size=1)
 a=runif(1,min=0,max=1)
 if((a<0.5)){
 if(i==1){return(c(vec[1]-step,vec[2:12]))}
 else if (i==12){return(c(vec[1:11],vec[12]-step))} 
 else
 {
 return(c(vec[1:(i-1)],vec[i]-step,vec[(i+1):12]))
 }
 }
 else if(vec[i]<9) 
 { 
 if(i==1){return(c(vec[1]+step,vec[2:12]))}
 else if (i==12){return(c(vec[1:11],vec[12]+step))}
 else 
 {
 return(c(vec[1:(i-1)],vec[i]+step,vec[(i+1):12]))
 }
 
 }
 
 }

 #Cruce entre dos cadenas de soluciones
 crossover<-function(r1,r2){
 i=sample(12,size=1) 
 if(i==1){return(c(r1[1],r2[2:12]))}
 else if (i==12){return(c(r1[1:11],r2[12]))}
 else 
 {return(c(r1[1:(i-1)],r2[i:12]))}
 }
 
 #Creación de la población inicial
 A=matrix(nrow=popsize,ncol=12)
 
 for(i in (1:popsize)){
 A[i,]=rbind(sample(9,size=12,replace=TRUE))
 }

 #Ganadores por generación
 topelite=as.integer(elite*popsize)
 
 #Loop principal
 for( i in (1:maxiter)){
 score=sapply(1:nrow(A),function(p){costf(A[p,])})
 B=cbind(A,score)
 Bdata<-data.frame(B)
 
 Bdata<-Bdata[with(Bdata,order(score)),]
 A=as.matrix(Bdata[1:topelite,1:12])
 while (length(A[,1])<popsize){
 a=runif(1,max=1,min=0)
 if(a<mutprod){
 c3=sample(topelite, size=1)
 A=rbind(A,mutacion(A[nrow=c3,1:12]))
 }
 else
 {
 c1=sample(topelite, size=1)
 c2=sample(topelite, size=1)
 A=rbind(A,crossover(A[nrow=c1,1:12],A[nrow=c2,1:12]))
 }
 } 
 
 }
 print(A[1,1:12])
 
 }

Se  compara el resultado y  se obtenemos lo siguiente.

#Método de Algoritmo Genético

 sol5
 6 5 4 5 6 7 4 5 5 7 5 5

costo.tiempo(sol5)
 2898

impri.horarios(sol5)
[1] "Seymour" "BOS" "13:40-15:37 138" "12:31-14:02 234"
[1] "Franny" "DAL" "10:30-14:57 290" "12:31-14:02 234"
[1] "Zooey" "CAK" "13:40-15:38 137" "15:07-17:21 129"
[1] "Walt" "MIA" "11:28-14:40 248" "12:31-14:02 234"
[1] "Buddy" "ORD" "12:44-14:17 134" "15:07-17:21 129"
[1] "Les" "OMA" "12:18-14:56 172" "12:31-14:02 234"

Es muy curioso como los resultados de los algoritmo muestran comportamientos muy distintos en cuando a las elecciones de vuelos y precios. Claro que esto es dependiente de la función de costo definida. El costo por medio de la solución de algoritmos genéticos es mucho mejor que las otras dos. Es casi 1000 unidades menos que la mejor de las otras dos.

Espero esta entrada de una idea de lo que es en general la optimización y el tipo de algoritmos que existen. Por su puesto que hay muchos temas que ni siquiera menciono, pero espero puedan consultarlos y que lo visto en esta entrada les ayude a tener una visión global del tema.

Comentario: El algoritmo que actualmente se emplea para hacer Machine Learning con grandes volúmenes de datos se llama Gradiente Descendiente Estocástico, se ha verificado que es adecuado para la mayoría de algoritmo de clasificación y predicción de ML, lo cual lo ha vuelto popular y se ha implementado en varias librerías, para ver  una breve, pero buena explicación el método recomiendo el vídeo de 12 minutos de Andrew Ng: https://class.coursera.org/ml-003/lecture/107

Referencias.

Libros:

http://www.amazon.com/Optimization-Machine-Learning-Information-Processing/dp/026201646X

http://www.princeton.edu/~sbubeck/Bubeck14.pdf

http://www.amazon.es/Stochastic-Optimization-Scientific-Computation-Schneider/dp/3540345590

Centros de investigación:

http://www.cs.manchester.ac.uk/mlo/

http://research.microsoft.com/en-us/groups/mlo/

Notas, artículos y cursos:

http://stoprog.org/SPTutorial/SPTutorial.php

http://www.people.fas.harvard.edu/~plam/teaching/methods/mcmc/mcmc_print.pdf

http://online.stanford.edu/course/convex-optimization-winter-2014

http://www.econ.uiuc.edu/~roger/research/conopt/coptr.pdf

Clasificación binaria….Naive Bayes

Un poco sobre análisis de textos

Lo principal de esta entrada es dar un ejemplo de como clasificar correos, esto lo hago siguiendo el ejemplo de libro Machine Learning for Hackers [1]. La técnica es conocida como Naive Bayes, pero creo conveniente comentar un poco de la librería de R para text mining.

La librería tm   permite realizar text mining en R project. El text mining se aboca sobre la conversión textos a datos que puedan ser analizados, ya sea por herramientas estadísticas, o por técnicas de procesamiento natural del lenguaje (NLP, en ingles).

La librería es bastante conocida en la red y existe mucha información, en varios blog se muestran ejemplos interesantes. En particular se mostraban análisis de textos de Twitter, un buen ejemplo es el de Yanchang Zhaog, este se puede encontrar en sus notas data mining con R. Lo desafortunado es que la API de twitter fue actualizada y ya no se puede replicar el ejemplo con el código de Yanchang.

Para ejemplificar el funcionamiento de tm tomo los títulos de las primeras planas de un periódico desde el día 1 al 23 de Marzo.

Los pasos son los siguientes:

  1. Extraemos los HTML de cada portada desde la página del periódico.
  2. Tomo el texto de la página web o HTML.
  3. Le quito a los textos la mayoría de palabras que forman parte de los títulos y secciones del periódico.
  4. Proceso con tm los archivos.
  5. Construyo una nubes de palabras.

Podríamos hacer otro tipo de análisis  después del punto 4, pero implica hablar de otros conceptos, como cluster, PCA  o redes. Observación: No pongo todo el código, pero prácticamente solo falta el proceso del paso 2.

#Paso 0
#Cargamos las librerías 
library(tm)
library(wordcloud)
library(stringi)

#Paso 1
#Extraemos los datos
direc<-""http://www.eluniversal.com.mx/hemeroteca/edicion_impresa_201503"
comp_dir<-".html"
dias<-c("01", "02", "03", "04", "05", "06", "07", "08", "09", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", "20", "21", "22", "23")
#Paso 2
#Extraemos los datos
edi_impresa<-character(length(dias))

for (i in 1:length(dias)){
     edi_impresa[i]=htmlToText(stri_paste(direc,dias[i],comp_dir))
    }
#Paso 3 y 4
#Construimos el Corpus
Doc<-Corpus(VectorSource(edi_impresa))

#Visualizamos Doc
Doc[[1]]

#Procesamos Doc con la librería tm

Doc1<-tm_map(Doc, function(x) stri_replace_all_regex(as.character(x), "<.+?>", " "))
Doc2 <- tm_map(Doc1, function(x) stri_replace_all_fixed(x, "\t", " "))
Doc3 <- tm_map(Doc2, PlainTextDocument)
Doc4 <- tm_map(Doc3, stripWhitespace)
Doc5 <- tm_map(Doc4, removeWords, stopwords("spanish"))
Doc6 <- tm_map(Doc5, removePunctuation)
Doc7<-tm_map(Doc6,removeNumbers)

#Hacemos una lista con todas las palabras de títulos o secciones del periódico
L<-#LO OMITO
Doc8 <- tm_map(Doc7, removeWords, L)
Doc9 <- tm_map(Doc8,content_transformer(tolower))

#Contruimos nuestra matriz de TDM
DocsTDM <- TermDocumentMatrix(Docs14)

Lo que hice hasta esta parte del código es la construcción de la matriz de términos de los documentos. Con ello puedo revisar las frecuencias de las palabras. Como observación, debido a que no hice bien la lista de palabras a remover aparecerán algunas que se refieren a títulos y secciones del periódico.

#Revisamos las frecuencias de las palabras
Grupo1<-rowSums(as.matrix(DocsTDM))
Grupo1<-subset(Grupo1,Grupo1>=25)

#Graficamos nuestras palabras
library(ggplot2)
qplot(names(Grupo1), Grupo1, geom="bar", xlab="Terms") + coord_flip()
barplot(Grupo1, las=2)

 Palabras_más_frecuentes

Por último hago la nube de palabras.

#Nube de palabras
Grupo1<-subset(Grupo1,Grupo1>=10)
m<-as.matrix(Grupo1)
word_freqs = sort(rowSums(m), decreasing=TRUE)
dm = data.frame(word=names(word_freqs), freq=word_freqs)
#Nube de palabras
wordcloud(dm$word, dm$freq, random.order=FALSE, colors=brewer.pal(8, "Dark2"))

Imagen_datos_nube

Bueno, la idea sería procesar correctamente los datos y presentar la nube de palabras con mayor frecuencia de apariciones en el texto o con algún criterio que de deseé. Otro ejemplo sería revisar cuales palabras persistieron por el paso de los días en el mes de Marzo. También se puede hacer otro tipo de análisis, por ejemplo buscar cluster  o tópicos en los datos, lo cual resulta más interesante.

Cómo observación, la librería tm cuenta con aproximadamente 54 funciones y cabe mencionar que el objeto principal es la construcción y procesamiento de Corpus, que es como un concentrado de todos los textos  que se analizan. Para procesar el contenido del Corpus se usa la función tm_map, la cual funciona análogamente a las funciones apply, lapply, sapplycon la diferencia de que se aplica solo a los Corpus. Para mayor detalle se puede revisar el manual de la librería [].

Si se tiene interés en aprender sobre procesamiento de textos siguiendo las metodologías de NLP, se tienen muy buenos cursos de la plataforma Coursera. Recomiendo el curso impartido por la Universidad de Stanford.

¿Web Scraping con R?

El procesar textos ha sido una área de investigación desde hace tiempo, ahora se tienen librerías y metodologías estándar, pero donde abunda mucho información para analizar y era complicado extraerla es en la web. Ha este tipo de programas se les llama web scraping.

Mario Annau liberó una librería en R llamada tm.plugin.webminig y por otro lado Hadley Wickham liberó la librería rvest. Las dos librerías se toman como modelo librerías de Python, las cuales son más robustas y se han mejorado con el tiempo.

No hago algún ejemplo, pero se puede revisar la documentación y replicar el ejemplo que comparten cada uno sobre el uso de su librería. Lo interesante es que las dos tratan de simplificar pasos que tm no hace por si sola, principalmente por la interacción con textos HTML.

Como recomendación ( y eso en cada entrada de esta categoría) creo importante revisar el repositorio de Hahley Wickham, el cual a contribuido a la comunidad de R project con librerías muy populares y de mucho utilidad, basta mencionar ggplot2 y reshape. Se ha convertido en el RockStar de R para muchos.

Clasificación de mensajes-detectando Spam

 

En la entrad sobre Máquina de Soporte Vectorial (SVM) hice un ejemplo comparativo para clasificar correos. Replico el código en esta entrada solo enfocándome en Naive Bayes.

#Se instala la librería kernlab para usar sus datos 
library(kernlab)
#Se extraen los datos
data(spam)
#Se genera una muestra de 100 índices
index=sample(1:nrow(spam),100)
#Se carga la librería e1071 para usar la función naiveBayes
library(e1071)
model=naiveBayes(type~.,data=spam[-index,])
model=naiveBayes(type~.,data=spam[-index,])
pred=predict(mol,spam[index,])
table(pred,spam[index,58])
 
pred   nonspam   spam
nonspam  38      0
 spam     29     33

Los datos que se usan en el código anterior son datos ya procesados, los cuales ya cuentan con la etiqueta o variable categórica de ser Spam o no. Lo recomendable es  cargar los datos y explorar la información, en este caso solo hago uso de los algoritmos Naive Bayes que se encuentra en la librería e1071.

Hago otro ejemplo con los datos de la librería kernlab, estos corresponden a 8993 personas de las cuales se toman 14 datos y uso el algoritmo de Naive Bayes para evaluar la clasificación de 1000 de estas personas por sexo.

#Se carga las librerías
library(kernlab)
library(e1071)
#Se extraen los datos
data("income")
#Se toman 1000 personas de los datos
index=sample(1:nrow(spam),1000)
#Se genera el modelo Naive Bayes de la librería e1071
modelo=naiveBayes(SEX~.,data=income[-index,])
#Se estima la predicción
pred=predict(modelo,income[index,])
#Se presentan los datos en modo de una tabla
table(pred,income[index,2])
 
pred  M   F
 M   234 162
 F   220 384

En este último ejemplo, al igual que el primero se hace uso de datos previamente procesados. Lo cual facilita la aplicación del algoritmo.

Las tablas de los dos ejemplos muestra dos categoría, los errores del algoritmo se pueden considerar como la cantidad de Hombres de la fila M de Male que son considerados como mujeres en la columna F de Famale, es decir; 384. Equivalente las 220 mujeres que son clasificadas como hombres en la segunda fila.

Estos errores son conocidos como falsos negativos y falsos positivos, el objetivo de las técnicas de clasificación es tener el mínimo de este tipo de errores.

El siguiente ejemplo consiste  también en clasificar correos, donde un conjunto son identificados como Spam y otro conjunto se sabe que no es spam. Estos son extraídos desde SpamAssasin, la dirección es http://spamassassin.apache.org/publiccorpus/ . Pero recomiendo bajar todos los datos y código libro Machine Learning for Hackers  dese el repositorio de Jhon M. White en Github, para replicar de manera completa el ejemplo.

A diferencia de los primero dos ejemplos de clasificación, en este caso los datos o textos se van a procesar para identificar las palabras que aparecen con mayor frecuencia en los Spam y las palabras de los que no son Spam.

Se hará uso de 3 funciones, la primera para extraer el mensaje del correo, la segunda para convertir el mensaje en una matriz de términos del documento, con la cual se revisa la frecuencias de las palabras y se aplica una tercera función para clasificar cualquier correo con respecto a las palabras que se encuentran con mayor frecuencia en los correos identificados como Spam.

Naive Bayes

Antes doy una breve explicación  de la idea a tras de la clasificación. Suponiendo que se tienen dos clase, hombres y mujeres; y una lista de personas, por ejemplo: Guadalupe, Ernesto, Sofía, Andrea, Guadalupe, Sofía, Luis, Miguel, Guadalupe y Jaime

Nombre     Sexo
Guadalupe- Mujer
Ernesto  - Hombre
Sofía    - Mujer
Andrea   - Mujer
Guadalupe- Hombre
Sofía    - Mujer
Luis     - Hombre
Miguel   - Hombre
Guadalupe- Mujer
Jaime    - Hombre

De la lista anterior al considerar dos categorías, hombres y mujeres, se puede tomar una persona, ejemplo Guadalupe,  pero uno puede preguntarse ¿cuál es la probabilidad de que sea mujer dado que la persona se llama Guadalupe?

P(Mujer|Guadalupe)=[P(Guadalupe|Mujer)P(Mujer)]/P(Guadalupe)

La ecuación anterior expresa P( | ) como la probabilidad condicional y la regla de Bayes. Con los datos del ejemplo anterior  se puede calcular el valor de P(Mujer|Guadalupe)=2/3.

Análogamente se puede calcular P(Hombre|Guadalupe)=1/3

Entonces  lo que se tiene es que tomando a una persona que se llama Guadalupe es más probable que sea mujer. Pero se pueden considerar más atributos, como edad, peso, color de ojos, etc. Así  se puede uno volver a preguntar, ¿cuál es la probabilidad de que sea Mujer dado que se llama la persona Guadalupe?, pero claro ahora se pueden considerar los nuevos atributos en el calculo de la probabilidad condicional.

Lo que se modificará en el calculo es la cantidad de factores, ejemplo:

P(Guadalupe|Mujer)=P(Peso|Mujer)*P(Color de ojos|Mujer)*P(Edad|Mujer)

En el ejemplo de los correos se tiene una lista de palabras que se encuentran en aquellos que son Spam, entonces el total de atributos será dado por el total de palabras. La pregunta que se hace es: ¿el correo es Spam dado que tienen estas palabras?

Lo que  no menciono en el calculo anterior es que los atributos  se consideran independientes condicionalmente. Los aspectos teóricos de manera breve se pueden consultar en la parte 4 de las notas de Andrew Ng o en las notas de Andrew Moore, con más detalles se pueden encontrar en la referencias [2,3,4].

Espero que la explicación pese a ser corta ayude a entender lo que se va hacer con los correos que se procesan.

Correos..¿spam o no?

Como antes mencioné, tomo un conjunto de correos identificados como Spam y otro que se tienen identificado que no lo es. Se procesan para formar una matriz de términos para cada conjunto de correos. Con la matriz se asignarán densidades a las palabras  y por último se construye una función para clasificar.

Los primeros 500 correos se consideran como el conjunto de entrenamiento para el proceso de clasificación.

#Cargamos los datos
library(tm)
library(ggplo2)
spam.path <- file.path("data", "spam")
easyham.path <- file.path("data", "easy_ham")

#Funciones que usaremos
######Primera función

get.msg <- function(path)
{
 con <- file(path, open = "rt")
 text <- readLines(con)
 msg <- text[seq(which(text =="")[1]+1, length(text), 1)]
 close(con)
 return(paste(msg, collapse = "\n"))
}

#####Segunda función

get.tdm2 <- function(doc.vec)
{
 library(stringi)
 doc.corpus<-Corpus(VectorSource(doc.vec))
 cor.doc<-tm_map(doc.corpus, function(x) stri_replace_all_regex(as.character(x), "<.+?>", " "))
 cor.doc <- tm_map(cor.doc, PlainTextDocument)
 cor.doc <- tm_map(cor.doc,content_transformer(tolower))
 cor.doc <- tm_map(cor.doc,removeWords,stopwords("english"))
 cor.doc<- tm_map(cor.doc, removePunctuation)
 cor.doc<- tm_map(cor.doc, removeNumbers)
 cor.doc<-tm_map(cor.doc, stripWhitespace)
 doc.dtm<- TermDocumentMatrix(cor.doc)
 return(doc.dtm)
}

######Tercer función
classify.email <- function(path, training.df, prior = 0.5, c = 1e-6)
{
 #Extrae el mensaje
 #Convierte los datos en una matrz de terminos del documento

 msg <- get.msg(path)
 msg.tdm <- get.tdm(msg)
 msg.freq <- rowSums(as.matrix(msg.tdm))
 msg.match <- intersect(names(msg.freq), training.df$términos)
 if(length(msg.match) < 1)
 {
 return(prior * c ^ (length(msg.freq)))
 }
 else
 {
 match.probs <- training.df$ocurrencia[match(msg.match, training.df$términos)]
 return(prior * prod(match.probs) * c ^ (length(msg.freq) - length(msg.match)))
 }
}

######Procesar los los Spam

spam.docs <- dir(spam.path)
spam.docs <- spam.docs[which(spam.docs != "cmds")]
all.spam <- sapply(spam.docs,function(p) get.msg(file.path(spam.path, p)))
spam.tdm <- get.tdm2(all.spam)

######Procesamos 500 correos que no son spam

easyham.docs<-dir(easyham.path)
easyhamo.docs500=easyham.docs[1:500]
easyham500<-sapply(easyham.docs500,function(p)get.msg(file.path(easyham.path,p))
easyham500.tdm<-get.tdm2(easyham500)

En el código anterior se obtiene las matrices de términos y se puede ver los valores de esta en R, solo colocando el nombre:

#Revisamos los valores de los corpus

spam.tdm

<<TermDocumentMatrix (terms: 23109, documents: 500)>>
Non-/sparse entries: 70715/11483785
Sparsity : 99%
Maximal term length: 298
Weighting : term frequency (tf)

#Corpus para correos que no son Spam
easyham500.tdm

<<TermDocumentMatrix (terms: 13118, documents: 500)>>
Non-/sparse entries: 50247/6508753
Sparsity : 99%
Maximal term length: 245
Weighting : term frequency (tf)

Con las matrices de términos se construye una matriz con las densidades de las palabras.

#Matriz de densidad

spam.matrix <- as.matrix(spam.tdm)
spam.counts <- rowSums(spam.matrix)
spam.df <- data.frame(cbind(names(spam.counts),as.numeric(spam.counts)),stringsAsFactors = FALSE)
names(spam.df) <- c("términos", "frecuencia")
spam.df$frecuencia <- as.numeric(spam.df$frecuencia)
spam.occurrence <- sapply(1:nrow(spam.matrix),
 function(i)
 {
 length(which(spam.matrix[i, ] > 0)) / ncol(spam.matrix)
 })
spam.density <- spam.df$frecuencia / sum(spam.df$frecuencia)

#Aplicamos el mismo proceso para easyham500.tdm

#Términos frecuentes en spam.df

head(spam.df[with(spam.df, order(-ocurrencia)),],n=15)
      términos frecuencia density ocurrencia
6019    email    837    0.007295962 0.574
15199   please   459    0.004001011 0.520
3539    click    350    0.003050880 0.454
11722   list     419    0.003652339 0.442
21468   will     843    0.007348262 0.442
7317    free     651    0.005674637 0.420
9990  information 364   0.003172915 0.374
13837   now      329    0.002867827 0.352
2899    can      518    0.004515302 0.344
7687    get      425    0.003704640 0.334
14323   one      376    0.003277517 0.322
13511   new      336    0.002928845 0.300
16452   receive  327    0.002850394 0.298
19437   time     316    0.002754509 0.296
12527   message  245    0.002135616 0.286

#Términos frecuentes en easy.df

 head(easy.df[with(easy.df, order(-ocurrencia)),],n=15)
       términos frecuencia density ocurrencia
5052    group   232      0.003446125 0.388
12391   use     272      0.004040284 0.380
12979   wrote   238      0.003535249 0.380
1577    can     348      0.005169187 0.368
6982    list    249      0.003698642 0.368
8258    one     358      0.005317727 0.338
6527    just    273      0.004055138 0.326
8120    now     231      0.003431271 0.324
4812    get     230      0.003416417 0.282
6924    like    232      0.003446125 0.282
3653    email   188      0.002792549 0.276
11337   subject 162      0.002406346 0.270
11854   time    188      0.002792549 0.258
12834   will    318      0.004723567 0.254
6112   information 162   0.002406346 0.232

 De los datos anteriores se observa que comparten los dos tipos de correos algunas palabras, pero se aprecia que las ocurrencias en los Spam tienen valores más altos.

Lo que sigue es comparar las palabras contra una muestra de correos que no son Spam, pero por el texto que contiene son más complicados para detectarse como no-Spam.

Usamos la función classify.email se prueba contra el conjunto de entrenamiento que son las matrices spam.df y easy.df

#Clasificación
hardham.path <- file.path("data", "hard_ham")
hardham.docs <- dir(hardham.path)
hardham.docs <- hardham.docs[which(hardham.docs != "cmds")]

hardham.spamtest <- sapply(hardham.docs, function(p) classify.email(file.path(hardham.path, p), training.df = spam.df))
 
hardham.hamtest <- sapply(hardham.docs,function(p) classify.email(file.path(hardham.path, p), training.df = easyham.df))

hardham.res <- ifelse(hardham.spamtest > hardham.hamtest,TRUE,FALSE)

summary(hardham.res)
   Mode FALSE TRUE NA's 
logical  242   7   0

Lo que se observa es que se clasifican 7 email como Spam, cuando no lo son. Es posible que se deba algún error en el procesamiento del TDM de los Spam, ya que el resultado obtenido es demasiado bueno.

Para cerrar esta entrada, hago lo que sería un paso de prueba del clasificador. Se inicia con un conjunto de entrenamiento y se prueba como responde para clasificar  2500 correos que no son Spam y que son fáciles de identificar (easyham).

#Función clasificador de Spam

spam.classifier <- function(path)
{
 pr.spam <- classify.email(path, spam.df)
 pr.ham <- classify.email(path, easy.df)
 return(c(pr.spam, pr.ham, ifelse(pr.spam > pr.ham, 1, 0)))
}
#Cargamos los 2500 correos

easyham.docs <- dir(easyham.path)
easyham.docs <- easyham.docs[which(easyham.docs != "cmds")]

#Los clasificamos

easyham.class <- suppressWarnings(lapply(easyham.docs,function(p){spam.classifier(file.path(easyham.path, p))}))

#Contruimos una matriz con los datos
easyham.matrix <- do.call(rbind, easyham.class)
easyham.final <- cbind(easyham.matrix, "EASYHAM")

class.matrix <- rbind(easyham.final)
class.df <- data.frame(class.matrix, stringsAsFactors = FALSE)

names(class.df) <- c("Pr.SPAM" ,"Pr.HAM", "Class", "Type")
class.df$Pr.SPAM <- as.numeric(class.df$Pr.SPAM)
class.df$Pr.HAM <- as.numeric(class.df$Pr.HAM)
class.df$Class <- as.logical(as.numeric(class.df$Class))
class.df$Type <- as.factor(class.df$Type)

#Hacemos la gráfica con la información de los correos clasificados

ggplot(class.df, aes(x = log(Pr.HAM), log(Pr.SPAM))) +
 geom_point(aes(shape = Type, alpha = 0.5)) +
 stat_abline(yintercept = 0, slope = 1) +
 scale_shape_manual(values = c("EASYHAM" =1),name = "Email Type") +
 scale_alpha(guide = "none") +
 xlab("log[Pr(HAM)]") +
 ylab("log[Pr(SPAM)]") +
 theme_bw() +
 theme(axis.text.x = element_blank(), axis.text.y = element_blank())

#Hacemos una tabla con las densidades 

easyham.col <- get.results(subset(class.df, Type == "EASYHAM")$Class)
class.res <- rbind(easyham.col)
colnames(class.res) <- c("NOT SPAM", "SPAM")
print(class.res)

          NOT SPAM SPAM
easyham.col 0.974 0.026

Se observamos que el clasificador identifica aproximadamente el 97.5% de los correos que no son Spam, lo cual es correcto y solo tienen un error (falsos positivos) de 2.5%. Este ejercicio se puede repetir con el resto de conjuntos de correos que se encuentran en los datos del texto Machine Learning for Hackers.

Grafica_Ham_vs_Spam

Espero que los ejemplos den una idea aproximada de lo que se hace con este algoritmo y como se usa para construir un clasificador, en el texto se tienen más detalles, pero en general traté de poner en esta entrada las ideas básicas.

También se debe de observar que la clasificación no se reduce solo a tener dos clases (Spam y no Spam), como mencioné en otras entradas existen varias técnicas de clasificación multi-clase y en muchas situaciones es probable que se tengan más de dos clases para analizar.

Referencias:

1.-http://shop.oreilly.com/product/0636920018483.do

2.-http://shop.oreilly.com/product/0636920018483.do

3.-http://www.amazon.com/Artificial-Intelligence-Modern-Approach-Edition/dp/0136042597

4.-http://statweb.stanford.edu/~tibs/ElemStatLearn/

5.-http://www.autonlab.org/tutorials/naive02.pdf

6.-http://cs229.stanford.edu/

Sobre sistemas de recomendación, ejemplo sencillo.

Sobre los sistemas de recomendación

Los sistemas de recomendación son quizás una de esas cosas con las cuales todos tenemos contacto actualmente, van desde sistemas sofisticados y predictivos, hasta otros menos rebuscados que solo nos hacen una lista fija de sugerencia. Los ejemplo conocidos son youtube, spotify, los sistemas de facebook, amazon, Netflix,etc.

Pese a que no es un tema fijo o algo que aparezca como tema dentro de las referencias importantes de Machine Learning, existen algunas referencias que hablan sobre el tema [2] y en particular un laboratorio reconocido por su investigación en sistemas de recomendación es GroupLens de la Universidad de Minnesota.

Este tipo de sistemas se hicieron famosos gracias a la competencia de Netflix,  en la cual un equipo logró cumplir con el reto de minimizar un indicador estadístico para asegurar que se obtenían mejores sugerencias [3]. Existen mucho artículos respecto a este concurso y sobre el sistema que se desarrollo, pero lo que se debe de resaltar es que el equipo ganador no hizo solo de un algoritmo de Machine Learning, uso varios y la parte más importante fue el pre-procesamiento de los datos [4].

Sobre el ejemplo

El siguiente ejemplo de sistema de recomendación lo presento como se discute el texto Machine Learning for Hackers [1]. En el ejemplo se hace uso de la técnica k-Nearest Neighbor(Knn) para definir las recomendaciones, es un método de clasificación no paramétrico. Antes de explicar un poco sobre Knn hago un ejemplo sencillo sobre series de tiempo y principalmente sobre la estimación de la “regresión de k-vecinos cercanos“.

#Regresión Vecinos Cercanos
Mortalidad<-scan("data/cmort.dat")
t=1:lenght(Mortalidad)
#Usamos la función supsmu
plot(t,Mortalidad,col="4",main="Mortalidad por semanas",xlab="Semanas",ylab="Número de muertes")
lines(supsmu(t,Mortalidad,span=0.5),col="2")
lines(supsmu(t,Mortalidad,span=0.01),col="4")

Ejemplo_Regresión_vecinos_cercanos

En la gráfica se observan dos curvas cercanas a los datos que están en color verde.Observando con detenimiento la curva en rojo se parece mucho a la línea de tendencia de los datos, pero no es tal cual una recta. Por otro lado la curva en azul, parece más cercana a una descripción del comportamiento “real” de los datos. En ambos casos lo que se hace es ir calculando una regresión para cierta cantidad de vecinos cercanos a cada punto, la línea roja se calcula para una cantidad mayor de vecinos en cada punto y la azul para una cantidad menor, por ello la línea roja no es una recta del todo y la azul tiene pequeños picos y a la vista nos parece mejor para describir los datos.

La gráfica anterior dice que cuando se calcula la regresión o el mejor ajuste lineal para una cantidad mayor de vecinos se tiende a seguir la tendencia de los datos. Cuando se calcula para una cantidad menor de vecinos es probable que la curva siga de cerca el comportamiento de los datos. Con este ejemplo solo deseo mostrar que Knn es hasta cierto punto útil para hacer exploración de la información y que si bien este caso es una regresión, la idea en general es la misma de considerar los vecinos cercanos. Esto de que “parezca” muy parecido a los datos tienen un costo  en el sentido estadístico en caso de considerar dichos algoritmos como candidatos para implementación.

Una imagen burda del algoritmo Knn es pensar en que si se está localizado en uno de los datos en color verde  se traza un circulo centrado en ese punto y se cuenta la cantidad k vecinos para luego hacer la regresión. Lo que se debe de resaltar es que trazar un circulo implica hablar de una distancia para comparar qué puntos están dentro o no del circulo, igual que en el análisis de Cluster se requiere la noción de distancia en el algoritmo de Knn.

Hago un ejemplo solo para mostrar como funciona la librería class que contiene la función knn que permite hacer el cálculo de k-vecinos cercanos. Tomo un subconjunto de entrenamiento y uno de prueba, estimo el modelo para el conjunto de datos de entrenamiento y comparo el resultado con los datos de prueba. Hago una comparación de un método lineal paramétrico de clasificación ( regresión logística) con respecto al método Knn.

Los datos empleados se pueden descargar desde github.

#Comparación de métodos
data.path<-file.path("data","example_data.csv")
datos<-read.csv(data.path)
n<-1:nrow(datos)
set.seed(1)

#Tomamos el 50% de los datos
library(class)
indices<-sort(sample(n,50))
training.x<-datos[indices,1:2]
training.y<-datos[indices,3]
test.x<-datos[-indices,1:2]
test.y<-datos[-indices,3]
#Aplicamos el método Knn
prediccion<-knn(training.x,test.x,training.y,k=5)
sum(prediccion!=test.y)
#Aplicamos un método lineal logit
modelo.logit<-glm(Label~ X+Y,data=datos[indices,])
prediccion2<-as.matrix(predict(modelo.logit,newdata=df[-indices,])>0)
sum(prediccion2!=test.y)

Se observa que en el ejemplo el método Knn da como predicción un 88% de coincidencias y el método lineal solo el 60%. Así que para situaciones no lineales se puede usar con mejor eficiencia el método Knn. Cuando digo predicción, me refiero a cuantos de los puntos de prueba el método clasificó correctamente y cuantos no.

Ejemplo de sistema de recomendación

Para el ejemplo se toma un conjunto de datos respecto a los paquetes instalados en R project. La intención de detectar de una lista de paquetes instalados cuando un usuario va  instalar otro paquete sabiendo cuales tiene instalados. Fue parte de una competencia en Kaggle. Lo interesante es ver que con una muestra pequeña se pudo aprender del comportamiento de los usuarios caracterizados por los paquetes instalados.

#Revisión de paquetes instalados
data.path<-file.path("data","isntallations.csv")
instalaciones<-read.csv(data.path)
#Revisamos las cabeceras
head(instalaciones)
#Construimos una matriz de datos
library(reshape)
matriz.de.paquetes<-cast(instalaciones,User~Package,value='Installed')
#Modificamos las filas 
row.names(matriz.de.paquetes) <- matriz.de.paquetes[, 1]
matriz.de.paquetes <- matriz.de.paquetes[, -1]
#Matriz de similitud
matriz.similitud <- cor(matriz.de.paquetes)
#Definimos una distancia
distancia <- -log((matriz.similitud / 2) + 0.5)

#Calculamos los vecinos cercanos
k.nearest.neighbors <- function(i, distancia, k = 25)
{
 return(order(distancia[i, ])[2:(k + 1)])
}

#Función para estimar la probabilidad
prob.instalacion <- function(user, package, matriz.de.paquetes, distancia, k = 25)
{
 neighbors <- k.nearest.neighbors(package, distancia, k = k)
 
 return(mean(sapply(neighbors, function (neighbor) {matriz.de.paquetes[user, neighbor]})))
}

#Función para elegir los paquetes más probables de isntalar
paquetes.mas.prob <- function(user,matriz.de.paquetes, distancia, k = 25)
{
 return(order(sapply(1:ncol(matriz.de.paquetes),
 function (package)
 {
 prob.instalacion(user,package,matriz.de.paquetes,distancia,k = k)
 }),
 decreasing = TRUE))
}

user <- 1

listing <- paquetes.mas.prob(user, matriz.de.paquetes, distancia)

colnames(matriz.de.paquetes)[listing[1:15]]
 [1] "adegenet" "AIGIS" "ConvergenceConcepts" "corcounts" "DBI" "DSpat" "ecodist" 
 [8] "eiPack" "envelope" "fBasics" "FGN" "FinTS" "fma" "fMultivar" 
[15] "FNN" 

Explico en lo general los pasos del código anterior. Primero  se cargan los datos que se encuentran en formato “csv”, luego se transforman a una matriz lo cual se hace con la función cast de la librería reshape. Una vez transformados, se calcula su matriz de correlaciones, pero la observación es que el método Knn no funciona con la matriz de correlaciones, sino con la matriz de distancias o similitudes.

El paso clave es usar el logaritmo y definir una distancia mediante la “normalización” de los valores obtenidos en la matriz de correlaciones, este paso es crucial para poder aplicar el algoritmo. Espero se pueda entender este paso, lo que se hace es considerar que las correlaciones tienen valores entre -1 y 1, suponiendo que no se alcanza el valor -1, entonce el intervalo de posibles valores tienen longitud 2. Si dividimos cada valor de la correlación sobre 2 y después se le sumamos 0.5, obligamos a que los valores estén entre 0 y 1. El logaritmo , de hecho el negativo del logaritmo con valores entre 0 y 1 indicarán que cuando los valores se aproximen a cero indica que son más distantes los paquetes de los usuarios  y que cuando se aproximen a 1 son muy cercanos.

Por último define una funciones que permite estimar la probabilidad mediante el cálculo de Knn. Entonces se obtiene una matriz la cual permite conocer cuales son los paquetes con mayor probabilidad para poder recomendarlos al usuario. La última línea muestra los 15 paquetes con mayor probabilidad de recomendar al usuario 1.

Entonces construir un sistema de recomendaciones es prácticamente tener una versión de este tipo de algoritmos o de este tipo de técnicas en el sistema, pero como lo mencioné al inicio de la entrada no existe un solo algoritmo útil para estos sistemas ya que se puede recurrir a varios para determinar cual funciona de manera adecuada. Ejemplo, si los datos muestran un comportamiento lineal es preferible usar la regresión logística que knn ya que clasificará mejor los datos.

Puede parecer confuso lo que hago en el ejemplo, pero en resumen es lo siguiente:

  1. Se cargaron los datos.
  2. Se construye una matriz de correlaciones con la información de los paquetes instalados en cada usuario.
  3. Se define una matriz de distancia por medio de la matriz de correlaciones.
  4. Se aplica el método Knn, el cual estimará los paquetes más cercanos a cada usuario.

Faltaría diseñar una interfaz agradable para mostrar los resultados.

Referencias:

1.-http://shop.oreilly.com/product/0636920018483.do

2.-http://www.mmds.org/

3.-http://www.wired.com/2009/09/bellkors-pragmatic-chaos-wins-1-million-netflix-prize/

4.-http://www.netflixprize.com/assets/ProgressPrize2008_BigChaos.pdf

Análisis de Cluster, un ejemplo sencillo.

El análisis de agrupamiento (cluster), es sencillo de explicar y de implementar en R project. Lo delicado es procesar la información y la definición de la métrica que consideramos apropiada para nuestros datos o experimento.

Hago un ejemplo sencillo antes de explicar qué se hace con esta técnica.

#Ejemplo de Cluster
data(iris)
#Generamos 5 grupos en nuestros datos
modelo5<-kmeans(iris[1:2],5)
plot(iris[1:2],col=modelo5$cluster,main="Ejemplo de 5 cluster",xlab="longitud de la Sápal",ylab="Ancho de la Sápal")
#Generamos 3 grupos en nuestros datos
modelo3<-kmeans(iris[1:4],3)
plot(iris[1:4],col=modelo3$cluster,main="Ejemplo de 3 cluster")

 Ejemplo_5_cluster

Ejemplo_3_cluster

Los dos gráficos muestran de colores distintos un conjunto de puntos, por lo que se ve en el código, usando una función kmeans para definir modelo5,modelo3. Luego indico al gráfico que tomara los colores partiendo de un valor $cluster. Esto podríamos deducirlo viendo solo el  código, la pregunta es cómo se define qué puntos pintar de cierto color. Si uno hace el “scatter plot” de los mismos datos se observa una “regla”.

#Scatter Plot 
data(iris)
plot(iris[1:4],col=c(2,3),main="Scatter Plot de parejas de vectores de datos")

Ejemplo_Scatter_Plot_Datos Iris

Si uno compara los 3 gráficos, la regla para determinar el color de los puntos tienen que ver con un modo de agruparlos, de tener un modo de verlos cuan cercanos o similares son. Eso se nota principalmente entre los dos últimos gráficos.

La idea de “similar”, puede ser pensada como un modo de comparar cosas, si se piensa en tres personas  donde dos estudiaron en la misma universidad, tienen casi la misma edad y alturas y tienen diferencia de peso de solo un par de kilos, son “casi iguales” con respecto al tercer individuo, que no estudió en la misma universidad, tiene más años, tiene sobrepeso y mide más que los otros dos. Entonces podemos pensar que los dos primeros son “similares” o al compararlos “son más cercanos o parecidos”.

La noción de “cercanía” está asociada con la noción de distancia, podemos definir que dos cafeterías tienen una distancia pequeña si las separa una calle, pensando que su unidad de medida son el número de calles entre ellas. Entonces la noción a través de todo esto es la de “métrica“, cómo medir similitudes.

En geometría, la métrica es la distancia entre puntos la cual tienen ciertas propiedades; medir del punto 1 al punto 2 es igual que iniciar midiendo desde el punto 2 hasta el punto 1 (simetría); la distancia entre los puntos es por lo menos cero (positiva) y  medir la distancia entre dos puntos es menor o igual que la suma de las distancia de los puntos a un punto intermedio (desigualdad del tringualo) y la última, si la distancia entre dos puntos es cero los puntos no pueden ser distintos (identidad de los indiscernibles).

Entonces la idea de la técnica es medir la distancia de cada pareja de puntos de los datos. Esto no necesariamente implica que la distancia sea para números o vectores, puede ser definida por las características de los datos o sobre cadenas de texto o mensajes, características de una persona,etc.

Aclarado el código anterior, lo que se hizo con el comando kmeans fue aplicar un algoritmo que permite separar entre los grupos de datos, tomando como parte inicial el  eligir la cantidad de grupos que se desean, después se calcula el centroide con respecto a todos los datos y se hace un proceso de separación e iteracciones para determinar los cúmulos, la técnica se llama k-Means Clustering. Para entender en detalles el algoritmo se puede revisar las notas de Andrew Ng o de Andrew Moore, las notas de éste último son acompañadas de varias  gráficas lo cual ayuda a visualizar el algoritmo de un modo sencillo.

Un ejemplo puede ser pensando en una matriz de datos, donde cada fila representa un individuo y cada columnas es una característica la cual puede tomar el valor de -1,0 o 1. Entonces lo que se hace es realizar el producto de la matriz y su transpuesta para obtener la correlación entre cada uno de los individuos y lo que se usa como métrica es la distancia euclidiana. Con ellos se puede tener una matriz con información sobre las distancia entre nuestros individuos.

#Ejemplos 
x.matriz<-matrix(sample(c(-1,0,1),24,replace=TRUE),nrow=4,ncol=6)
#Etiquetamos las columnas y las filas
row.names(ex.matriz)<-c('A','B','C','D')
colnames(ex.matriz)<-c('P1','P2','P3','P4','P5','P6')
#Multiplicamos por la transpuesta
ex.mult<-ex.matriz%*%t(ex.matriz)

#Calculamos la distancia
ex.dis<-dist(ex.mult)

Ahora usamos el comando cmdscale para calcular el escalamiento multidimensional, obtendremos como resultado los puntos que son cercanos o próximos, según la métrica.

#Cluster
ex.mds<-cmdscale(ex.dis)
plot(ex.mds,type='n')
text(ex.mds,c('A','B','C','D'))

Ejemplo_Escalado_Multidimensional

Se concluye del ejemplo anterior que la estrategia es; procesar nuestra información, definir o usar una métrica, construir la matriz de similaridad y aplicar el análisis de escalamiento multidimensional o un algoritmo de agrupamiento.

Polarización de un Senado

Muestro un ejemplo del texto Machine Learning for Hackers [1]. El ejemplo usa los datos del senado de EEUU, se encuentran en la plataforma http://www.voteview.com y en caso de tener una falla en la descargar se puede bajar del repositorio de GitHub en https://github.com/johnmyleswhite/ML_for_Hackers.

La idea es la siguiente, revisar el comportamiento de las votaciones del Senado y revisar el comportamiento del periodo 110 al 111 del congreso.

Lo que se hará es procesar la información, asignar ciertas variables categóricas que solo tendrán valores -1,0 y 1 para cada una de las características de los senadores analizados. Esto permite definir una métrica sencilla, como la del último ejemplo y basta considerar el productos de matrices para formar la matriz de similaridad.

#Análisis de votos 
library(foreign)
library(ggplot2)
#Tomamos la ruta de los archivos
data.dir 6, 0,no.pres[,i])
         no.pres[,i] 0 & no.pres[,i] < 4, 1, no.pres[,i])
         no.pres[,i] 1, -1, no.pres[,i])
 }
 
 return(as.matrix(no.pres[,10:ncol(no.pres)]))
}

#Aplicamos la función a cada entrada de rollcall.data
rollcall.simple

Con lo anterior pueden construir las gráficas que para apreciar la evolución del senado, donde la hipótesis inicial es que el congresos se polarizó en el 111 congreso.

Primero se revisa el comportamiento solo de las votaciones del 111 congreso y después se compara con respecto los congresos previos, partiendo desde el 101.

MDS_Votacion_Senado111_US

MDS_Votacion_Senado del 101 al 111_congreso_USEn resumen la última gráfica muestra que el Senado en distintos periodos se puede dividir en dos grupos, se puede decir que se polariza entre Republicanos y Demócratas. Los pocos senadores considerados como “independientes” parece, al revisar la última gráfica; que siempre votaran como los Demócratas.

La moraleja es que bajo una métrica sencilla se pueden extraer ciertas propiedades del comportamiento o dinámica del los datos. Sería interesante hacer un ejercicio similar para  el Senado de otros países.

Resumiendo, la técnica o análisis de agrupamientos (cluster) requiere como primer paso explorar nuestra información y procesarla si es necesario para poder definir una métrica que tenga coherencia con el tipo de información y datos.

Después se debe de calcular la matriz de distancias y por último aplicar escalamiento multidimensional o algún algoritmos para detectar aglomeraciones.

Como paso intermedio, tienen sentido plantearse una pregunta respecto a los posibles grupos que se busca encontrar en los datos y como siempre para cerrar el análisis, hacer usos de gráficos que permitan visualizar el resultado.

Ha esta entrada le falta mucho muchos ejemplos, que espero más adelante compartir otros. Respecto al tema falta mostrar varios detalles; tipos de métricas para tipos de datos, algoritmos de aglomeración para tipos de análisis y la creación de modelos basados en esta técnica.

Los aspectos teóricos se encontrarán en la categoría “Sobre Machine Learning”. Una disculpa para la gente de matemáticas por no haber puesto la definición formal de espacio métrico, métrica y mencionar algunas de sus propiedades, espero posteriormente hacer algunas entradas con un enfoque formal.

Actualización 19-05-2015

En la primera versión de esta entrada faltó mencionar muchos detalles sobre la técnica de K-medias, pero además no mencioné nada sobre otro algoritmo Hirarchical Clustering o Cluster jerárquico.

Comparto un ejemplo de esta técnica y en resumen lo que se hace en este algoritmo es construir un árbol que ayude aglomerar los datos en conjuntos los cuales al terminar el algoritmo muestra las aglomeraciones obtenidas.

El costo computacional es mucho mayor que los métodos k-medias, pero en general para muestras de datos relativamente grandes no resulta tan costoso implementarlo y además se han desarrollado algoritmos para reducir el costo en tiempo y procesamiento.

Los datos los tomo de la página del libro Elements of Stadistical Learning [2], en el código indico como descargarlo directamente.

Respecto al origen de los datos, cabe resaltar que corresponde al análisis de tumores cancerígenos, se analiza su valor genético a 64 tumores y por medio de las técnicas de Cluster se trata de detectar cuantos cánceres corresponde a qué tipo de tumor.

Algunas precisiones. Se descargaron los nombres de los tumores en un archivo txt, a que los nombres de las columnas no vienen por default en los datos. La matriz de datos está formada por 64 columnas y 6830 filas.

#Ejemplo de Cluster con datos de tumores 
#Ubicado en el directorio donde estan los nombres de los genes se procede a cargarlos en R
#Se carga como un vector
Nombres=read.table("nombres_tumores.txt")
#Visualizamos los datos
head(Nombres)
# V1
# CNS
# CNS
# CNS
# RENAL
# BREAST
# CNS

#Cargamos los datos de los genes

datos=scan("http://statweb.stanford.edu/~tibs/ElemStatLearn/datasets/nci.data")

#Visualizamos los datos
head(M)
#Tranformamos los datos en una matriz de 64 renglones por 6830 columnas
#Después invertimos la matrix
M=matrix(datos,nrow=64)
M2=t(M)

#Agremagamos los nombres a las columnas de M2

colnames(M2)=as.vector(t(Nombres))
 
head(M2)
#Así se visualiza la tabla de datos
# CNS CNS CNS RENAL BREAST CNS CNS BREAST NSCLC NSCLC RENAL RENAL RENAL RENAL RENAL RENAL RENAL BREAST NSCLC RENAL UNKNOWN
#[1,] 0.300 0.679961 0.940 2.800000e-01 0.485 0.310 -0.830 -0.190 0.460 0.760 0.270 -0.450 -0.030 0.710 -0.360 -0.210 -0.500 -1.060 0.150 -0.290 -0.200
#[2,] 1.180 1.289961 -0.040 -3.100000e-01 -0.465 -0.030 0.000 -0.870 0.000 1.490 0.630 -0.060 -1.120 0.000 -1.420 -1.950 -0.520 -2.190 -0.450 0.000 0.740
#[3,] 0.550 0.169961 -0.170 6.800000e-01 0.395 -0.100 0.130 -0.450 1.150 0.280 -0.360 0.150 -0.050 0.160 -0.030 -0.700 -0.660 -0.130 -0.320 0.050 0.080
#[4,] 1.140 0.379961 -0.040 -8.100000e-01 0.905 -0.460 -1.630 0.080 -1.400 0.100 -1.040 -0.610 0.000 -0.770 -2.280 -1.650 -2.610 0.000 -1.610 0.730 0.760
#[5,] -0.265 0.464961 -0.605 6.250000e-01 0.200 -0.205 0.075 0.005 -0.005 -0.525 0.015 -0.395 -0.285 0.045 0.135 -0.075 0.225 -0.485 -0.095 0.385 -0.105
#[6,] -0.070 0.579961 0.000 -1.387779e-17 -0.005 -0.540 -0.360 0.350 -0.700 0.360 -0.040 0.150 -0.250 -0.160 -0.320 0.060 -0.050 -0.430 -0.080 0.390 -0.080

#Aplicamos k-Means con 3 cluster a
KM3=kmeans(M2,3)

#Generamos una gráfica para visualizar los componentes

library(cluster)
clusplot(M2,KM3$clusteres,shape=TRUE,color=TRUE,labels=3)

#Aplicamos Hierarchical Clustering
H1<-hclust(dist(M),"average")
#Dendrograma
plot(H1,col=2,frame.plot=TRUE)

 Las gráficas que obtengo del código  son las  siguientes:

Dendrograma

Cluster_Genes

La primera gráfica es el modo visual en que se comportan los clusters en el método jerárquico y la segunda es una herramienta visual para observar el comportamiento del método k-medias.

Una observación entre los dos métodos es, k-medias pide como parámetro conocer cuantos clusters se buscamos en el algoritmo y el jerárquico no, de cierto modo la aglomeración por el jerárquico es cercana a lo que sería un método no supervisado.

Para el  algoritmo de k-medias resulta importante el análisis de la elección de la cantidad de clusters en los datos.

Existen varias técnicas para tratar de elegir el mejor valor de k, en R project, existen varias librerías entre las cuales están: “fpc” y “mclust”. Para los datos de este ejemplo probé con mclust, pero le resultaba muy costoso así que cancelé el proceso.

Pero otro modo muy sencillo y que no requiere ninguna librería es mediante el análisis de los cuadrados de las distancias de cada  dato aglomerado con respecto a la cantidad de clusters que se calculan.

#Elección de K para k-medias
wss <- (nrow(M2)-1)*sum(apply(M2,2,var))
 for (i in 2:10) wss[i] <- sum(kmeans(M2,centers=i)$withinss)
plot(1:10, wss, type="b", xlab="Númbero de Clusters",ylab="Suma de cuadrados dentro de cada cluster")

 La gráfica que se obtiene es la siguiente:

K-Medias

Esta última técnicas es sencilla, pero el problema es que no siempre indica de modo claro como elegir la cantidad de clusters de un modo óptimo. Por lo cual es bueno recurrir a otro tipo de algoritmo y compara nuestra selección y al final validar con alguna muestra de datos.

En un trabajo de investigación de Robert Tibshirani-G. Walther y Trevor Hastie[5] plantean un método para obtener un óptimo del número de clusters en el algoritmo k-meadias, sobre este reporte de investigación se basan varias librerías en R project para estimar el valor óptimo de k.

Espero que este complemente de una idea en general del tipo de análisis y de las técnicas básicas de como usar este tipo algoritmo. Cabe mencionar que las técnicas de Cluster están relacionadas con la Mezcla de Modelos y el algoritmo EM, lo cual se puede usar para la construcción de modelos. En otra entrada explico el uso y relación de esto.

Referencias:

1.-http://shop.oreilly.com/product/0636920018483.do

2.-http://statweb.stanford.edu/~tibs/ElemStatLearn/

3.-http://www.amazon.com/Cluster-Analysis-Brian-S-Everitt/dp/0470749911

4.-http://www.amazon.com/Machine-Learning-Probabilistic-Perspective-Computation/dp/0262018020

5.-http://web.stanford.edu/~hastie/Papers/gap.pdf

Regresión no lineal, Cross-Validation y Regularization.

En esta entrada sigo la postura del texto Machine Learning for Hackers[1], la cual es explicar los tres temas juntos para mostrar como se conectan. Se pueden revisar por separado de manera amplia y desde una perspectiva estadística en la referencia [2].

En breve, se puede mencionar que Cross-Validation y Regularization funcionan bien usándolos para determinar un buen modelo o determinar cuando se está “sobre estimando” o “sobre ajustado” (overfitting).

La regresión no-lineal, como su nombre lo indica; es además de una generalización de la regresión lineal, también otra opción para construir un modelo. Al ser más general, no es extraño que se pueda caer en la idea de pensar que cuanto más complejo el modelo es mejor, lo cual no es necesariamente cierto.

Algo que suele ocurrir es que uno busca por medio de un modelo lineal estudiar un fenómenos no lineales, o en otro caso implementar alguna técnica que funciona bien con fenómenos no lineales en fenómenos que pueden ser “bien” estudiados con modelos lineales. Por lo cual es importante explorar la información y tratar de construir más de un modelo, lineal o no-lineal y hacer uso de herramientas para determinar cual es un mejor modelo.

Sobre la idea de complejidad de un modelo, se pude hacer un análisis bastante serio al respecto, pero uno se acercaría al concepto de complejidad  y este puede tener varios puntos de vista los cuales dejo sin detalle en esta entrada, basta mencionar que la complejidad tienen varias aristas y puntos de vista, que van desde temas de matemática puras hasta temas de sistemas complejos y teoría de computación.

Hago un  primer  un ejemplo de ajuste lineal y no-lineal.

#Ejemplo no lineal
x<-seq(-10,10,by=0.01)
y<-1-x^2+rnorm(length(x),0,7)
ggplot(data.frame(X=x,Y=y),aes(x=X,y=Y))+geom_point(col=2)+geom_smooth(size=3,method='lm',se=FALSE)
ggplot(data.frame(X=x,Y=y),aes(x=X,y=Y))+geom_point(col=2)+geom_smooth(size=3,se=FALSE)

Ajuste_LinealAjuste_No-Lineal

La gráfica muestra que al calcular de manera lineal la regresión, la recta estimada corta los datos y al calcular la regresión no-lineal se estima una curva, la cual parece un mejor ajuste. Pero si pensamos en transformar los datos de entrada, elevándolos al cuadrado, se obtiene otro tipo de comportamiento.

#Ejemplo de ajusto con datos transformados
x.cuadrados<-x^2
ggplot(data.frame(X.Cuadrados=x.cuadrados,Y=y),aes(x=X.Cuadrados,y=Y))+geom_point(col=2)+geom_smooth(size=3,se=FALSE)+ ggtitle("Ajuste de datos transformados")+theme(plot.title = element_text(lineheight=.8, face="bold"))

 

Ajuste_No_lineal_con_datos_transformados

La gráfica anterior muestra casi una relación lineal con el ajuste no-lineal, se puede revisar los residuos o la r cuadrada para estimar cuanto de la variabilidad de Y es explicada con los ajustes. Lo hago en R  por medio del comando summary y el resultado es el siguiente:

#Revisamos las estimaciones
> summary(lm(y~x))$r.squared
[1] 0.0001486085
> summary(lm(y~x.cuadrada))$r.squared
[1] 0.9444159

En el ejemplo se observa como una transformación en los datos hace que se pase de un ajuste que explica prácticamente el 0%, a un ajuste que explica casi el 95% de la variabilidad de los datos de salida.

Del ejemplo se puede concluir, que una transformación  hace una gran diferencia y en este caso beneficia al modelo, pero lo que no se visualiza es el costo de transformar los datos, nada limita a pensar que se  puede considerar cualquier potencia de la información y revisar el ajuste. Con potencia me refiero a elevar a la tercera, cuarta o quinta potencia nuestros datos de entrada.

Antes de abordar el problema de determinar cuando un ajuste no-lineal es mejor que otro, hago otro ejemplo para mostrar los ajustes no-paramétricos. Para este ejemplo uso datos que pueden descargarse desde la página de Robert H.Shumway.

#Serie de tiempo de la mortalidad cardiovascular
mortalidad<-scan("cmort.dat")
t=1:length(mortalidad)
lm(mortalidad~t)
t2=t^2
t3=t^3
c=cos(2*pi*t/52)
s=sin(2*pi*t/52)
fit=lm(mortalidad~t+t2+t3+c+s)
#Graficamos el ajuste
plot(t,mortalidad,type="p",main="Datos de Mortalidad Cardiovascular",xlab="Tiempo",ylab="Cantidad de muertes",lwd="2")
lines(fit$fit,col="2")
abline(a=96.65,b=-0.0312,col="3")
lines(lowess(t,mortalidad),col="4")
lines(ksmooth(t,mortalidad,"normal",bandwidth=104),col="5")
lines(ksmooth(t,mortalidad,"normal",bandwidth=5),col="6")

Regresiones_lineales_y_nolineales

La gráfica del ejemplo anterior muestra varios ajustes, el primero es el lineal, el segundo el polinomial y el tercero y cuarto son métodos no paramétricos que pueden ser considerados como extensiones de la regresión polinomial (regresión no-lineal). Cada uno muestra cierto comportamiento, y puede servir para explorar la información y posteriormente pensar en el tipo de modelo que se puede construir. Con la gráfica anterior solo trato mostrar que las regresiones no se limitan solo a lo lineal y no-lineal, uno puede hacer uso de técnicas no-paramétrica para estimar otro tipo de regresiones. Para mayor detalles se puede revisar la presentación de Andrew W. Moore [3].

Lo que hago ahora es enfocarme a la regresión polinomial, construyo un ejemplo por medio del cual comparo el ajuste lineal y el polinomial.

#Ejemplo no lineal
x<-seq(0,1,by=0.01)
y<-sin(2*pi*x)+rnorm(length(x),0,0.1)
df<-data.frame(X=x,Y=y)
ggplot(df,aes(x=X,y=Y))+geom_point()+ggtitle("Relación no lineal")+theme(plot.title = element_text(lineheight=.8, face="bold"))

Relación_No_Lineal

Se observa en la gráfica que los datos parecen como si estuvieran en una montaña rusa, tienen partes altas y bajas. Pareciera que siguen una trayectoria curva y queda pensar en aplicar dos regresiones, una lineal y otra polinomial.

#Ajuste lineal y polinomial
summary(lm(Y~X,data=df))
Call:
lm(formula = Y ~ X, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.91351 -0.39942  0.02522  0.39833  0.98606 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.94573    0.09068   10.43   |t|)    
(Intercept)  -0.13568    0.04621  -2.936  0.00415 ** 
X            11.50556    0.40215  28.610  < 2e-16 ***
X2          -33.75407    0.93713 -36.018  < 2e-16 ***
X3           22.55399    0.61585  36.622  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1204 on 97 degrees of freedom
Multiple R-squared:  0.9726,    Adjusted R-squared:  0.9717 
F-statistic:  1147 on 3 and 97 DF,  p-value: < 2.2e-16

De los ajuste se  puede observar que el lineal explica aproximadamente el 59% y el polinomial con un polinomio de grado 3, explica aproximadamente el 97%. Pero nada limita a considerar polinomios de grado mayor, hago la estimación con un polinomio de grado 14 y  comparo como se comporta el ajuste, para esto usamos la función poly de R.

#Ajuste lineal y polinomial
summary(lm(Y~poly(X,degree=14),data=df))
Call:
lm(formula = Y ~ poly(X, degree = 14), data = df)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.243497 -0.066481 -0.001563  0.066558  0.240361 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)             0.004371   0.009707   0.450    0.654    
poly(X, degree = 14)1  -5.516374   0.097551 -56.549  < 2e-16 ***
poly(X, degree = 14)2   0.058760   0.097551   0.602    0.549    
poly(X, degree = 14)3   4.410330   0.097551  45.211  < 2e-16 ***
poly(X, degree = 14)4  -0.034386   0.097551  -0.352    0.725    
poly(X, degree = 14)5  -0.708767   0.097551  -7.266 1.58e-10 ***
poly(X, degree = 14)6  -0.005587   0.097551  -0.057    0.954    
poly(X, degree = 14)7   0.155736   0.097551   1.596    0.114    
poly(X, degree = 14)8   0.147905   0.097551   1.516    0.133    
poly(X, degree = 14)9   0.022788   0.097551   0.234    0.816    
poly(X, degree = 14)10  0.055761   0.097551   0.572    0.569    
poly(X, degree = 14)11 -0.087618   0.097551  -0.898    0.372    
poly(X, degree = 14)12  0.003417   0.097551   0.035    0.972    
poly(X, degree = 14)13  0.147284   0.097551   1.510    0.135    
poly(X, degree = 14)14 -0.075299   0.097551  -0.772    0.442    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.09755 on 86 degrees of freedom
Multiple R-squared:  0.984,     Adjusted R-squared:  0.9814 
F-statistic: 378.9 on 14 and 86 DF,  p-value: < 2.2e-16

Se observa que el ajuste explicado por el polinomio de grado 14 es del 98% de variabilidad en los datos, por cual puede hacer pensar que es mejor modelos, comparando los 3 que realicé.

Pero nada hace especial el polinomio de grado 14 con respecto a cualquier otro, bien podrías preferir uno de grado 20, 30 o 50.

Si reviso la R cuadrada de varios modelos polinomiales de grado mayor que 14 obtendré que explican del 98% en adelante y claro sin llegar  explicar el 100%. Pero el problema que surge tienen que ver con la siguiente pregunta,¿cómo elegir el mejor modelo de entre un conjunto de modelos?

En Machine Learning la idea no es encontrar un solo modelo, es construir varios, ya que más de una técnica puede funcionar para el mismo tipo de situación. Pero resulta crucial tener un modo o método para  elegir entre todos los modelos posibles.

Cross-Validation y Regularization

“All models are wrong, but  some models are userful”.

George Box

Este tema puedo ser consultado para más detalles y breve en las famosas nota de Andrew Ng, pero también de modo breve puede ser revisado en las notas de Andrew W. Moore. Recomiendo las notas del último, son más amigables ya que las de Andrew Ng requieren formación en matemáticas (probabilidad, álgebra lineal y cálculo de varias variables). Para ver estos dos temas de manera amplia en la referencia[2]

La idea de Cross-Validation es la siguiente:

  1. Elegir dos conjuntos, el de entrenamiento (training set) y el de prueba (test set).
  2.  Estimar el modelo en los datos de entrenamiento.
  3. Calcular MSE o RMSE de los datos de prueba.
  4. Comparar entre los Modelos con datos del paso 3.

En general se elige de todos los datos un 80-20 o 70-30,  es decir; 80% de los datos para entrenamiento y 20% para prueba. Pero se pueden tener algunos problemas con esto, el principal es que si el modelo requiere muchos datos por ser no-paramétrico es posible que la estrategia falle si no se cuenta con suficientes datos. Entonces se hace uso de una técnica conocida como LOOCV (Leave-one-out Cross-Validation), que en lugar de fijar el 20% para prueba solo no se considera uno de los datos, se calcula el modelo y se mide el error del dato no considerado y se hace ese proceso para toda la información, por último se calcula MSE o RMSE con la información de los errores.

Puede parecer que LOOCV cumple bien ante las situaciones donde la cantidad es importante, pero se puede usar otra alternativa, en lugar de hacer solo para una muestra de 20% fija, se divide la muestra en una cantidad finita de subconjuntos. Suponemos en K subconjuntos disjuntos, entonces se hace la estimación del modelo con cada uno de los subconjuntos tomados como el conjunto de prueba y se procede a calcular MSE o RMSE. Esta técnicas se llama K-Fold Cross Validation

La regularización está íntimamente relacionada con el concepto de complejidad en un modelo, se puede pensar que tener un modelo de un polinomio de grado 3 es más complejo que un polinomio de grado 2. Los fundamentos teóricos están relacionados con la noción de un “espacio normado“. No explico a que se refiere el último concepto, lo que interesa en esta entrada es entender cómo validar que un modelo es más complejo que otro.

Como la idea no es desarrollar la teoría, entonces lo que hago es usar una librería en R que realiza la estimación de los modelos lineales generalizados (GLM). Parte de sus algoritmos de estimación se realiza por regularización, así se puede usar uno de sus parámetros de estimación como medida de regularización.

Antes de hacer el ejemplo en R, lo que hago es mostrar como estimamos un modelo de regresión polinomial con respecto a unos datos y analizar el parámetro de regularización por medio de cross-validation para determinar cual es el mejor modelo. Considero que desea usar como modelo la regresión de un polinomio de grado 10, así se elige un conjunto de entrenamiento (training set) y uno de prueba (test set).

#Ejemplo del uso de regularización y cross-validation
x<-seq(0,1,by=0.01)
y<-sin(2*pi*x)+rnorm(length(x),0,0.1)
n<-lenght(x)
indices<-sort(sample(1:n,round(0.5*n))
training.x<-x[indices]
training.y<-y[indices]
test.x<-x[-indices]
test.y<-y[-indices]
df<-data.frame(X=X,Y=y)
training.df<-data.frame(X=training.x,Y=training.y)
test.df<-data.frame(X=test.x,Y=test.y)
#Usamos el conjunto de entrenamiento
glmnet.fit<-with(training.df,glmnet(poly(X,degree=10),Y)
lambdas<-glmnet.fit$lambda
modelo<-dta.frame()
for( lambda in lambdas){
modelo<-rbind(modelo,data.frame(Lambda=lambda,RMSE=rsme(test.y,with(test.df,predict(glmnet.fit,poly(X,degree=10),s=lambda)))))
}
#Visualizamos la relación entre Lambdas vs RMSE
ggplot(modelo,aes(X=Lambda,Y=RMSE))+geom_point()+geom_line()
#Elegimos Lambda que minimiza RMSE
best.lambda<-with(modelo, Lambda[which(RMSE==min(RMSE))])
glmnet.fit<-with(df,glmnet(poly(X,degree=10),Y))
#Visualizamos los coeficientes estimados del "mejor" modelos
coef(glmnet.fit,s=best.lambda

Del código anterior y del proceso que se realiza, se puede concluir que entre las dos herramientas Regularización y Cross-Validation se puede elegir el modelo que tienen menos “desviación estándar” en la distribución de los errores.

Se aprecia en la gráfica de abajo que el mínimo de la curva  ayuda a determinar el mejor modelo.

Lambda_vs_RMSE

 Caso de estudio

El siguiente ejemplo forma parte de los ejemplos presentados en el texto de Machine Learning for Hackers[1]. Lo que se hace en el texto es tomar datos de los 100 libros más populares de la editorial O’Reilly y tratar de definir partiendo de los textos de las contra portadas el posible éxito o alta venta.

Para desarrollar el ejemplo se usan las librerías tm y glmnet. Por medio de la primera se construye una matriz mediante los textos de la contra portada y ya teniendo la matriz se usa la regresión lineal general contra una variable para tratar de definir un rango de ventas partiendo de la lista de palabras que  contiene.

Mediante de la regularización (lambda) y la cross-validation, se trata de elegir el modelo que minimice RMSE con respecto al parámetro lambda del modelo.

#Análisis de contra-portadas
rango<-read.csv(data.file,stringsAsFactor=FALSE)
documentos<-data.frame(Texto=rango$Long.Desc.)
row.names(documentos)<-1:nrow(documentos)
library(tm)
corpus<-Corpus(DataframeSource(documentos))
corpus<-tm_map(corpus, content_transformer(tolower))
corpus<-tm_map(corpus,stripWhitespaces)
corpus<-tm_map(corpus,removeWords,stopwords('english'))
dtm<-DocumentTermMatrix(corpus)

rendimiento<-data.frame()

for(lambda in c(0.1,0.25,0.5,1,2,2.5))
{ 
  for ( i in 1:50){

    indices<-sample(1:100,80)
    training.x<-x[indices,]
    training.y<-y[indices]
    test.x<-x[-indices,]
    test.y<-y[-indices]
    glmnet.fit<-glmnet(training.x,training.y)
    prediccion.y<-predict(glm.fit,test.x,s=lambda)
    rmse<-sqrt(mean((prediccion.y-test.y)^2))
    rendimiento<-rbind(rendimiento,data.frame(Lambda=lambda,Interaccion=i,RMSE=rmse))
    }
}

El código anterior hace los siguiente:

  1. Se define el rango del parámetro de regularización (lambda).
  2. Se eligen subconjunto, 80% para entrenamiento y 20% para probar.
  3. Se estima la regresión y se compara contra el conjunto de prueba, con esta comparación se calcula RMSE para cierto lambda.

Análisis_Text_RMSE vs Lambda

Análisis_Text_RMSE vs Lambda2

 

La primera gráfica muestra que no se logra definir un buen modelo. En general RMSE vs Lambda nos indican que no se puede obtener mucho de los datos analizados, es decir; no resulta evidente como poder distinguir entre los 50 libros más vendidos con ese modelo.

Por lo cual se hace uso de una variable dummy para pasar de una regresión a una clasificación, por medio de la cual se puede responder cuando un libro es considerado uno de los 50 más vendidos solo partiendo de su contra portada. En este ejemplo no se trata solo de hacer el modelo, sino de verificar mediante regularización y cross-validation como es posible elegir un modelo con menor complejidad o uno mejor entre los posibles

Las observaciones y pasos a seguir son los siguientes; al usar la variable dummy y aplicar una regresión logística, se eligen subconjuntos de entrenamiento y prueba. Después de varias interacciones se puede analizar el comportamiento del parámetro lambda y  en lugar de RMSE  se usa la tasa de errores entre las predicciones y el subconjunto de prueba.

#Análisis de contra-portadas
y <- rep(c(1, 0), each = 50)
regularizacion.fit <- glmnet(x, y, family = 'binomial')
predict(regularizacion.fit, newx = x, s = 0.001)
ifelse(predict(regularizacion.fit, newx = x, s = 0.001) > 0, 1, 0)

performance <- data.frame()

for (i in 1:250)
{
 indices <- sample(1:100, 80)
 
 training.x <- x[indices, ]
 training.y <- y[indices]
 
 test.x <- x[-indices, ]
 test.y <- y[-indices]
 
 for (lambda in c(0.0001, 0.001, 0.0025, 0.005, 0.01, 0.025, 0.5, 0.1))
 {
 glm.fit <- glmnet(training.x, training.y, family = 'binomial')
 
 predicted.y <- ifelse(predict(glm.fit, test.x, s = lambda) > 0, 1, 0)
 
 tasa.error <- mean(predicted.y != test.y)

 performance <- rbind(performance,data.frame(Lambda = lambda,Iteration = i,TasaError = tasa.error))
 }
}

ggplot(performance, aes(x = Lambda, y = ErrorRate)) +
 stat_summary(fun.data = 'mean_cl_boot', geom = 'errorbar') +
 stat_summary(fun.data = 'mean_cl_boot', geom = 'point') +
 scale_x_log10()

 Análisis_Text_Tasa_Error vs log_ Lambda

En conclusión, se observa que al realizar el proceso con una regresión logística, se puede analizar el comportamiento de lambda y la tasa de error. Por lo tanto se puede clasificar los textos y elegir un modelo con la tasa de error correspondiente a lambda menor (como en la gráfica anterior).

¿Por qué el último modelo permite hacer una elección entre ellos y el anterior no?

El último modelo muestra constancia en la tasa de error, entonces se puede considerar que la menor Lambda hace considerar al modelo como el mejor.

En el modelo previo no se veía un patrón definido, elegir el que tiene menor RMSE o la menor Lambda no asegura que se este haciendo una buena elección.

Lo anterior puede parecer superficial, pero pensándolo un poco no es  una mala estrategia proceder con una variable dummy para modificar el modelo y tomar las tasas de error en lugar de RMSE.

Espero la entrada de una idea general de como se usa la regularización y la cross-validation, también como el modelo no-lineal puede beneficiarnos en algunos casos. En otra ocasión espero extender los ejemplo y hacer otro que puedan ayudar a clarificar los conceptos.

Algunas cosas que desearía fortalecer de esta entrada más adelante es resaltar con otro ejemplo como la regresión logística toma un rol de método de clasificación lineal, el rol de la regularización y la validación-cruzada o cross-validation para elegir modelos. Un par de buenos artículos sobre el tema se pueden ver en las referencias [2,5,6,7].

 Referencias:

1.-http://shop.oreilly.com/product/0636920018483.do

2.-http://statweb.stanford.edu/~tibs/ElemStatLearn/

3.-http://openclassroom.stanford.edu/MainFolder/DocumentPage.php?course=MachineLearning&doc=exercises/ex5/ex5.html

4.-http://www.amazon.com/Machine-Learning-Probabilistic-Perspective-Computation/dp/0262018020

5.-http://robjhyndman.com/hyndsight/crossvalidation/

6.-http://stats.stackexchange.com/

7.-http://projecteuclid.org/download/pdfview_1/euclid.ssu/1268143839