Ejemplo de Procesamiento en R y en Python

Previos.

Hace algunos meses en una competencia en Kaggle se presentó un código para un sistema de recomendación. En general para todo sistema de recomendación la parte más pesada suele ser el procesamiento, es quizás la clave de obtener buenas variables para tener un modelo más robusto. En fin, se pueden combinar muchas técnicas de todo tipo para crear un sistema de recomendación, en otro momento escribiré sobre lo poco que se al respecto.

El sistema de recomendación en el fondo terminaba siendo un problema de clasificación. Hubo diversos acercamiento para abordar el problema, todos con ciertos aspectos en común pero otros radicalmente distintos, les recomiendo leer la descripción de las soluciones las cuales las pueden encontrar aquí.

Lo interesante del código que tomé para esta entrada, me parece que es un buen código, por que es  legible, hacía operaciones de procesamiento estándar, un modelo “sencillo” y daba buen rendimiento ante la métrica de la competencia. El reto que pensé y el origen de esta entrada, fue replicar el procesamiento que se realizaba en R pero ahora en el entorno de Python.

En la migración del código respeté todos los pasos seguidos en el original, hice uso de la memoria de manera similar y trate de hacer la eliminación de tablas de manera igual. Esto para tratar de comparar lo mejor posible el uso de los dos entornos.

La lucha, R vs Python.

El debate entré los enamorados de un entorno y otro, suelen tener discusiones de todo tipo, en general no aportan nada y los argumentos a favor o en contra no tienen ni sentido ni valor alguno al momento de abordar o resolver un problema real.

Nunca el lenguaje será lo más importante si tienes buenas ideas, buenos algoritmos, etc. Hay gente que programa sus redes neuronales directo desde cuda con C o C++, sin pasar por tensorflow o pytorch. Así que esta entrada no es para apoyar un lenguaje u otro, es para comparar algunos aspectos:

  • Legibilidad del código.
  • Eficiencia en líneas de código.
  • Uso de la memoria.

Quizás hay otros criterios o detalles que se pueden comparar, pero creo que con esos basta.

Más detalles pueden ser considerados, pero creo que para limitar y mantener la entrada a un nivel intermedio los tres aspectos cubren lo general del código.

Código

Pego todo el código para su lectura. El código original en R fue escrito por Fabien Vavrand y la versión en Python es mía.

###########################################################################################################
#
# Kaggle Instacart competition
# Fabien Vavrand, June 2017
# Simple xgboost starter, score 0.3791 on LB
# Products selection is based on product by product binary classification, with a global threshold (0.21)
#
###########################################################################################################

Codigo en R

library(data.table)
library(dplyr)
library(tidyr)


# Load Data ---------------------------------------------------------------
path <- "../input"

aisles <- fread(file.path(path, "aisles.csv"))
departments <- fread(file.path(path, "departments.csv"))
orderp <- fread(file.path(path, "order_products__prior.csv"))
ordert <- fread(file.path(path, "order_products__train.csv"))
orders <- fread(file.path(path, "orders.csv"))
products <- fread(file.path(path, "products.csv"))


# Reshape data ------------------------------------------------------------
aisles$aisle <- as.factor(aisles$aisle)
departments$department <- as.factor(departments$department)
orders$eval_set <- as.factor(orders$eval_set)
products$product_name <- as.factor(products$product_name)

products <- products %>% 
 inner_join(aisles) %>% inner_join(departments) %>% 
 select(-aisle_id, -department_id)
rm(aisles, departments)

ordert$user_id <- orders$user_id[match(ordert$order_id, orders$order_id)]

orders_products <- orders %>% inner_join(orderp, by = "order_id")

rm(orderp)
gc()


# Products ----------------------------------------------------------------
prd <- orders_products %>%
 arrange(user_id, order_number, product_id) %>%
 group_by(user_id, product_id) %>%
 mutate(product_time = row_number()) %>%
 ungroup() %>%
 group_by(product_id) %>%
 summarise(
 prod_orders = n(),
 prod_reorders = sum(reordered),
 prod_first_orders = sum(product_time == 1),
 prod_second_orders = sum(product_time == 2)
 )

prd$prod_reorder_probability <- prd$prod_second_orders / prd$prod_first_orders
prd$prod_reorder_times <- 1 + prd$prod_reorders / prd$prod_first_orders
prd$prod_reorder_ratio <- prd$prod_reorders / prd$prod_orders

prd <- prd %>% select(-prod_reorders, -prod_first_orders, -prod_second_orders)

rm(products)
gc()

# Users -------------------------------------------------------------------
users <- orders %>%
 filter(eval_set == "prior") %>%
 group_by(user_id) %>%
 summarise(
 user_orders = max(order_number),
 user_period = sum(days_since_prior_order, na.rm = T),
 user_mean_days_since_prior = mean(days_since_prior_order, na.rm = T)
 )

us <- orders_products %>%
 group_by(user_id) %>%
 summarise(
 user_total_products = n(),
 user_reorder_ratio = sum(reordered == 1) / sum(order_number > 1),
 user_distinct_products = n_distinct(product_id)
 )

users <- users %>% inner_join(us)
users$user_average_basket <- users$user_total_products / users$user_orders

us <- orders %>%
 filter(eval_set != "prior") %>%
 select(user_id, order_id, eval_set,
 time_since_last_order = days_since_prior_order)

users <- users %>% inner_join(us)

rm(us)
gc()


# Database ----------------------------------------------------------------
data <- orders_products %>%
 group_by(user_id, product_id) %>% 
 summarise(
 up_orders = n(),
 up_first_order = min(order_number),
 up_last_order = max(order_number),
 up_average_cart_position = mean(add_to_cart_order))

rm(orders_products, orders)

data <- data %>% 
 inner_join(prd, by = "product_id") %>%
 inner_join(users, by = "user_id")

data$up_order_rate <- data$up_orders / data$user_orders
data$up_orders_since_last_order <- data$user_orders - data$up_last_order
data$up_order_rate_since_first_order <- data$up_orders / (data$user_orders - data$up_first_order + 1)

data <- data %>% 
 left_join(ordert %>% select(user_id, product_id, reordered), 
 by = c("user_id", "product_id"))

rm(ordert, prd, users)
gc()


# Train / Test datasets ---------------------------------------------------
train <- as.data.frame(data[data$eval_set == "train",])
train$eval_set <- NULL
train$user_id <- NULL
train$product_id <- NULL
train$order_id <- NULL
train$reordered[is.na(train$reordered)] <- 0

test <- as.data.frame(data[data$eval_set == "test",])
test$eval_set <- NULL
test$user_id <- NULL
test$reordered <- NULL

rm(data)
gc()


# Model -------------------------------------------------------------------
library(xgboost)

params <- list(
 "objective" = "reg:logistic",
 "eval_metric" = "logloss",
 "eta" = 0.1,
 "max_depth" = 6,
 "min_child_weight" = 10,
 "gamma" = 0.70,
 "subsample" = 0.76,
 "colsample_bytree" = 0.95,
 "alpha" = 2e-05,
 "lambda" = 10
)

subtrain <- train %>% sample_frac(0.1)
X <- xgb.DMatrix(as.matrix(subtrain %>% select(-reordered)), label = subtrain$reordered)
model <- xgboost(data = X, params = params, nrounds = 80)

importance <- xgb.importance(colnames(X), model = model)
xgb.ggplot.importance(importance)

rm(X, importance, subtrain)
gc()


# Apply model -------------------------------------------------------------
X <- xgb.DMatrix(as.matrix(test %>% select(-order_id, -product_id)))
test$reordered <- predict(model, X)

test$reordered <- (test$reordered > 0.21) * 1

submission <- test %>%
 filter(reordered == 1) %>%
 group_by(order_id) %>%
 summarise(
 products = paste(product_id, collapse = " ")
 )

missing <- data.frame(
 order_id = unique(test$order_id[!test$order_id %in% submission$order_id]),
 products = "None"
)

submission <- submission %>% bind_rows(missing) %>% arrange(order_id)
write.csv(submission, file = "submit.csv", row.names = F)

¿Largo y confuso?…calma, explicaré lo que me resulta importante.

Pero ahora veamos el mismo código completo en Python.

###########################################################################################################
#
# Kaggle Instacart competition
# Similary to Fabien Vavrand's script , Ago 2017
# Simple xgboost starter, score 0.3791 on LB
# Products selection is based on product by product binary classification, with a global threshold (0.21)
# 
# Daniel Legorreta 
# 
###########################################################################################################

import numpy 
import pandas
from sklearn.preprocessing import LabelEncoder
import xgboost as xgb

# Load Data ---------------------------------------------------------------
path="../input/"

aisles=pd.read_csv(path+"aisles.csv")
departments=pd.read_csv(path+"departments.csv")
orderp=pd.read_csv(path+"order_products__prior.csv")
ordert= pd.read_csv(path+"order_products__train.csv")
orders= pd.read_csv(path+"orders.csv")
products= pd.read_csv(path+"products.csv")

# Reshape data ------------------------------------------------------------

Factor=LabelEncoder()
Factor.fit(aisles.aisle)
aisles['aisle']=Factor.transform(aisles.aisle)
Factor.fit(departments.department)
departments['department']=Factor.transform(departments.department)
Factor.fit(products.product_name)
products['product_name']=Factor.transform(products.product_name)

products=departments.join(aisles\
 .join(products.set_index("aisle_id"),how="inner",on='aisle_id')\
 .set_index("department_id"),how="inner",on='department_id')

del(products['aisle_id'])
del(products['department_id'])
del(aisles,departments)

ordert=pd.merge(ordert,orders[orders.eval_set=='train'][['order_id','user_id']],how='left',on='order_id')
orders_products=pd.merge(orders,orderp, how='inner',on = 'order_id')

del(orderp)

# Products ----------------------------------------------------------------

Aux_4=orders_products[['user_id','order_number','product_id','reordered']]\
 .assign(product_time=orders_products\
 .sort_values(['user_id','order_number','product_id'])\
 .groupby(['user_id','product_id'])\
 .cumcount() + 1)

prd1=Aux_4.groupby('product_id')\
 .apply(lambda x:pd.Series(dict(prod_orders=x.user_id.count(),prod_reorders=x.reordered\
 .sum(),prod_first_orders=(x.product_time == 1).sum(),prod_second_orders=(x.product_time == 2).sum())))

prd1.loc[:,'prod_reorder_probability']=prd1.prod_second_orders/prd1.prod_first_orders
prd1.loc[:,'prod_reorder_times']=1+prd1.prod_reorders/prd1.prod_first_orders
prd1.loc[:,'prod_reorder_ratio']=prd1.prod_reorders/prd1.prod_orders

prd=prd1.drop(['prod_reorders', 'prod_first_orders', 'prod_second_orders'],axis=1)

del(Aux_4,prd1)


# Users -------------------------------------------------------------------


users=orders[orders['eval_set'] == "prior"]\
 .groupby('user_id')\
 .agg({'order_number':'max','days_since_prior_order':['sum','mean']})

users.columns = ["_".join(x) for x in users.columns.ravel()]
users.columns=["user_orders","user_period","user_mean_days_since_prior"]

us=orders_products[['user_id','reordered','order_number','product_id']]\
 .groupby('user_id')\
 .apply(lambda x:pd.Series(dict(
 user_total_products=np.size(x.product_id),user_reorder_ratio = np.divide((x.reordered == 1).sum(),(x.order_number > 1).sum()),user_distinct_products =np.size(x.product_id.unique()))))

users1=pd.merge(users,us, left_index=True, right_index=True,how='inner')
users1.loc[:,'user_average_basket']=users1.user_total_products / users1.user_orders

del(us)

us=orders[orders.eval_set != "prior"][['user_id', 'order_id', 'eval_set','days_since_prior_order']]\
 .rename(index=str, columns={"user_id": "user_id", "order_id": "order_id","eval_set":"eval_set","time_since_last_order":"days_since_prior_order"}) 
users=pd.merge(users1,us,left_index=True,right_on="user_id")

del(us,users1)


# Database ----------------------------------------------------------------

data=orders_products[["user_id", "product_id","order_number","add_to_cart_order"]]\
 .groupby(["user_id", "product_id"])\
 .agg({'order_number':['min','max','size'],'add_to_cart_order':['mean']})

data.columns = ["_".join(x) for x in data.columns.ravel()]
data.reset_index(level=[0,1], inplace=True)
data.columns =["user_id","product_id","up_first_order","up_last_order","up_orders","up_average_cart_position"]

prd.reset_index(level=0, inplace=True)

data=users.join(prd\
 .join(data.set_index("product_id"),on='product_id',how="inner")\
 .set_index("user_id"),on="user_id",how="inner")

data.loc[:,"up_order_rate"]=data.up_orders / data.user_orders
data.loc[:,"up_orders_since_last_order"]=data.user_orders - data.up_last_order
data.loc[:,"up_order_rate_since_first_order"]=data.up_orders / (data.user_orders - data.up_first_order + 1)

data=pd.merge(data,ordert[["user_id", "product_id", "reordered"]],how='left',on=["user_id", "product_id"])

del(ordert,prd,users)


# Train / Test datasets ---------------------------------------------------

train=data[data.eval_set == "train"]
train=train.drop(labels=['eval_set','user_id','product_id','order_id'], axis=1)
train.reordered=train.reordered.fillna(0)

test=data[data.eval_set == "test"]
test=test.drop(labels=['eval_set','user_id','reordered'], axis=1)

del(data)

# Model -------------------------------------------------------------------
import xgboost as xgb


subtrain=train.sample(frac=0.5)

param={'objective':'reg:logistic','eval_metric':'logloss','eta':0.1,'max_depth':6,'min_child_weight':10,
'gamma':0.7,'subsample':0.76,'colsample_bytree':0.95,'alpha':2e-05,'lambda':10}

X_train =xgb.DMatrix(subtrain.drop( "reordered", axis=1), label = subtrain.loc[:,"reordered"])


num_round = 80
model = xgb.train(param, X_train, num_round)

X_test =xgb.DMatrix(test.drop(labels=['order_id','product_id'],axis=1))

test.loc[:,'reordered']=model.predict(X_test)
test.loc[:,'reordered']=np.where(test.reordered>.21,1,0)


test.loc[:,'product_id']=test.product_id.apply(lambda x: str(x))
Submission=test[test.reordered==1][['order_id','product_id']].groupby('order_id')['product_id']\
 .agg(lambda x: ' '.join(x)).reset_index()

missing=pd.DataFrame()
missing.loc[:,'order_id']=np.unique(test.order_id[~test.order_id.isin(Submission.order_id)])
missing.loc[:,'product_id']=None

test_salida = pd.concat([Submission, missing], axis=0)
test_salida.columns=['order_id','products']

#Out writing
test_salida.to_csv("submit/submit4.csv", index=False)

¿Abrumado?…calma, revisemos los detalles que creo que son importantes.

Lo primero a comentar es que los entornos son sencillos, solo requiere unas cuantas bibliotecas para procesar los datos:

  • Por parte de R: data.table, dplyr, tidyr y xgboost.
  • Por parte de Python: numpy, pandas, sklearn y xgboost.

El que sean unas cuantas bibliotecas creo que me parece un buen ejemplo para comparar los dos entornos. Solo se hace uso de lo mínimo en los dos lenguajes, lo estándar para procesar y para construir un modelo de árboles con xgboost. Posiblemente algunos se confundan con algunas partes, si desean entender mejor todas las secciones podría explicar a detalle  pero creo innecesario para esta entrada.

Legibilidad.

Para se un poco justos, considerando que se tienen un nivel básico de R deberían resultar familiares el uso y funcionalidad de los “símbolos”:$,<-,%>%.

En Python debería resultar familiar la asignación con “=”, el uso de “.” y “\”.

Considerando que se tiene ese conocimiento mínimo, uno puede leer el código, pese a que no se conozcan muy bien las funciones y operaciones que se realizan, a mi parecer desde la parte de “#Reshape data” se empieza hacer un poco más difícil de seguir el código de Python.

Pero la situación se agrava cuando pasamos a leer parte del código donde se hace operaciones del tipo Filter+GruopBy+Agragation. Ejemplo en el siguiente extracto de código de la sección “# Users” en R:

users <- orders %>%
 filter(eval_set == "prior") %>%
 group_by(user_id) %>%
 summarise(
 user_orders = max(order_number),
 user_period = sum(days_since_prior_order, na.rm = T),
 user_mean_days_since_prior = mean(days_since_prior_order, na.rm = T)
 )

Para Python:

users=orders[orders['eval_set'] == "prior"]\
 .groupby('user_id')\
 .agg({'order_number':'max','days_since_prior_order':['sum','mean']})

users.columns = ["_".join(x) for x in users.columns.ravel()]
users.columns=["user_orders","user_period","user_mean_days_since_prior"]

Este tipo de operaciones entran en el marco de “split-apply-combine”, sobre este tipo de operaciones hace tiempo escribí unas entradas con ejemplos de como se utilizan y cual es la idea detrás, se pueden leer aquí. Para más detalles de este tipo de operaciones se puede leer los artículos [1,2].

Las operaciones suelen ser muy usadas en sql (SELECT+WHERE+GROUPBY), lo que observo es que pese a que los dos entornos pueden hacer lo mismo , en R se vuelve más legible este tipo de procesamientos y en Python resulta ser menos claras, necesitas tener en mente otro tipo de conceptos como estructuras de datos básicas, diccionarios,funciones lambda, funciones en Numpy, etc.

Si observamos con cuidado, se hacen las siguientes operaciones: Orden->Filtro->GroupBy->Agregación. Existen en pandas la función “filter”, pero no funcionaría similar al filtro condicional que necesitamos. La agrupación es similar en los dos entornos, pero cuando pasamos a crear las nuevas variables o las “agregaciones” resulta fácil entender el código de R pero la parte de Python termina de parecer un poco extraño.

Si bien las dos últimas líneas son para agregar un prefijo y renombrar las columnas, esto podría realizarse de otro modo, pero de igual forma resulta más “oscuro” que el código de R.

El siguiente fracción en R y Python creo muestra otro aspecto que muestra lo confuso que puede ser Python.

us <- orders_products %>%
 group_by(user_id) %>%
 summarise(
 user_total_products = n(),
 user_reorder_ratio = sum(reordered == 1) / sum(order_number > 1),
 user_distinct_products = n_distinct(product_id)
 )

En Python:

us=orders_products[['user_id','reordered','order_number','product_id']]\
 .groupby('user_id')\
 .apply(lambda x:pd.Series(dict(
 user_total_products=np.size(x.product_id),user_reorder_ratio = np.divide((x.reordered == 1).sum(),\
 (x.order_number > 1).sum()),\
 user_distinct_products =np.size(x.product_id.unique()))))

En mi opinión, es fácil pensar en elegir algunas columnas( filter), agrupar los datos con respecto alguna de las variables (groupby) y construir nuevas variables (summarize) que informen aspectos como la cantidad de objetos, la media, etc. Lo cual debería ser fácil traducir la idea a código.

En Python resulta complejo pensar el proceso, debido a que tienes que pensarlo con un nivel de “programación mayor” que el requerido en R. Tampoco es que sean necesarios aspectos sumamente avanzados de programación, pero pensando en la situación de trabajo estándar en R creo que te focalizas más en el proceso y en Python tienes que pensar un poco más en la programación, además de pensar en el proceso.

En varios pasos en Python tienes que pensar en aspectos delicados como los índices de las tablas y aplicar alguna operación sobre ellos. Esto es bueno y malo, puede frustrar en un inicio pero después con algo de práctica son útiles.

En resumen las operaciones son de un tipo estándar en SQL o el manejadores de bases de datos, el entorno de R resulta fácil y te libera para pensar en procesamiento y dejar un poco de lado la programación, en Python creo que te obliga a contar con una idea de programación más demandante que en R.

¿Cuanto de esto es culpa mía? 

Quizás toda, ya que yo escribí la versión en Python, pero confieso que me esforcé en tratar de escribirlo lo más simple y legible y respetar la estructura original del código en R.

Cantidad de Líneas de Código y Eficiencia.

Si copias y pegas el código, veras que aproximadamente es considerablemente más compacto el código en Python, por aproximadamente unas 30 líneas de código.

Podría parecer poco, pero si piensas en que tienes que hacer 10 códigos similares, estamos hablando de 300 líneas menos. Lo cual puede ser bastante ahorro de trabajo.

En buen parte el uso de “%>%” en R y de “\” en Python, simplifican mucho el trabajo, otro aspecto es que Pandas en si hace cierta programación que se asemeja a la programación funcional y el uso de “.” ahorra muchos pequeños pasos.

La eficiencia, pensándola en el uso de la memoria RAM y del tiempo en el procesamiento resultaba considerablemente mejor en Python. Cuando se cargan las tablas y se hace el procesamiento resultaba mucho mejor manejada la memoria, si observas en el código de Python trato de seguir el mismo uso de memoria para que fuera lo más similar posible.

En mi experiencia, Python me da mejor rendimiento que R cuando trabajo con tablas considerablemente grandes y cuando lanzo algún algoritmo suele ser aún más notorio el rendimiento. Esto es mi caso, quizás a otros les resulta mejor el entorno de R.

Un aspecto extra.

Los join, la secuencia de join o merge, pese a que se pueden realizar las mismas operaciones en los dos entornos me sorprendió que cuando hacer una sucesión de estas operaciones, el orden es el contrario en R que en Python. Ejemplo en la siguiente fracción de código:

products <- products %>% 
 inner_join(aisles) %>% inner_join(departments) %>% 
 select(-aisle_id, -department_id)

Para Python:

products=departments.join(aisles\
 .join(products.set_index("aisle_id"),how="inner",on='aisle_id')\
 .set_index("department_id"),how="inner",on='department_id')

El orden de operaciones es Productos -> aisles->departaments, al final es un inner join de las tres tablas. Cuando realizar esta cadena de operaciones haces uso de los indices de cada tabla, para poder combinarlas,

Se observa que el orden es al revés, podría parecer que fue mi culpa, pero al construir y comparar las tablas que se obtienen la secuencia de operaciones debían de seguir ese orden por la relación entre los índices.

Esto yo lo veo del siguiente modo, la cadena de operaciones en R van de afuera hacia dentro, por otro lado en Python para de adentro hacia afuera.

Puede que exagere, pero mi apreciación es que ciertas operaciones en Pandas son pensadas solo para un par de objetos, cuando pasar de esa cantidad se vuelve menos claro y poco intuitivo. Creo que en Pandas puede resultar quizás mentalmente complejo pensar en esas operaciones de manera “natural”.

Creo que en este código se repite lo mismo, el nivel de codificación requerida en Python resulta mayor que en R, que sea bueno o malo, no lo se, depende de nuestra formación y acercamiento ambos entornos.

Conclusión

El ejercicio de hacer el símil de un lenguaje a otro creo que siempre es bueno, ayuda a ver ciertos detalles nos parecen obvios. En lo personal el rendimiento del código en Python me sorprendió, resultaba muchísimo mejor que el código en R. Desconozco los detalles a bajo nivel como para saber si es por el tipo de operaciones entre tablas, el manejo de los índices o si terminar pagando con rendimiento el que sea más legible que se gana en R.

Espero te sirva el ejemplo para comparar el mismo tipo de operaciones y no está de más quizás hacer el ejercicio en un entorno de scala con spark.

Referencias:

  1. The Split-Apply-Combine Strategy for Data Analysis.
  2. Split-Apply-Combine en Pandas.
  3. Los datos se pueden descargar desde aquí.
  4. Los códigos se pueden descargar desde el repositorio dlegor.
Anuncios

¿Cuánto se puede saber desde los discursos?

Esto es un juego.

Esta entrada la estuve pensando en la última semana de Septiembre 2015, en especial debido a varios acontecimientos y fechas que se celebran en México. Ejemplo, el 1° de Septiembre se hace entrega del informe anual del gobierno por parte del presidente, el día 15 se celebra el día de la independencia de México, el día 27 de Septiembre del 2014 sucedieron unos hechos lamentables en el estado de Guerrera dónde 43 estudiantes de una Normal Rural fueron “desaparecido”. Pero más aún, con motivo de este último suceso se han realizado manifestaciones y posiblemente se vuelven desde este año actos que se realizarán año con año.

Lo que observé en Facebook, en twitter  y en la prensa, suele ser variado dependiendo del perfil tanto de las personas como de los periódicos. Pero algo que observé fue que las críticas se centraban en los discursos del presidente, sobré algo que decía, sobre lo dicho en la ONU con sede en New York mientras se llevaban acabo manifestaciones en México. Sobre el tipo de palabras, sobre los comentarios, etc.

Algunas frases fueron tomadas como los objetivos centrales de las críticas de dichos discursos. Entiendo que existen un equipo que redacta los discursos del presidente, pero también entiendo que existe un personal que revisa y que verifica el tipo de términos y explicaciones que se dirán. Evito imaginar que se tenga una estrategia para hacer crítica o freno a las aspiraciones políticas de un posible candidato a la presidencia, cuando aún faltan 3 años para la contienda electoral, la entrada la hago sin tener una postura política sobre las conclusiones e interpretaciones que se dieron a los discursos.

Pensé que no sería mala idea jugar con una muestra pequeña de discursos (48) y analizar cosas básicas para tratar de ver cómo ha cambiado la “similaridad” entre los discursos, ver cuáles han sido los tópicos a los largo de los meses y tratar de ver qué se puede saber de manera estadística desde los discursos. Pienso que son otros los capacitados para hacer un análisis de tiempo y espacio donde se dice tal o cual discurso y la relevancia o repercusión de lo dicho.

Esto no desentona en tipo de entradas, ya que no es para nada una manifestación política o una entrada donde trato de mostrar algo negativo del gobierno o donde trato de concluir o deducir algo sobre lo que sucede en el país. Es simplemente jugar con un puñado de datos y ver que se descubre.

Con toda mi ignorancia sobre como interpretar un discurso, me atrevo a suponer de un modo simplista que cada discursos político tienen algo de “localidad” y de “temporalidad”. Con localidad, pienso que depende de donde es emitido, de cuan tan importante es el lugar donde se emite y lo temporal me refiero al momento que rodea a dicho discurso, los acontecimiento cercanos, los sucesos sociales y políticos que acontecen en las fechas en las que se dice dicho discurso. Por otro lado, imagino que otro factor que afecta la “relevancia” de un discurso es gracias a los medios de comunicación, ya que bien un discurso puede ser parte de un evento donde un mandatario comunica algo o puede ser usado como parte de la información que usa y discute un medio de comunicación, el cual termina perneando la opinión publica.

Esto de cierto modo me permite hacer una análisis de textos, clasificar los discursos, detectar tópicos en los discursos, comparar por medio de medidas de similaridad como han cambiado o en qué meses se muestra mayor similaridad, comparar la muestra con respecto al informe y la participación en la ONU, etc. Estos aspectos pueden ser trabajados con una combinación de herramientas y algoritmos, lo cual pensé que sería divertido ver hasta donde se puede saber algo desde los discursos.

Sobre la muestra.

Tomé 48 discursos, 4 por mes desde Septiembre 2014 hasta Septiembre 2015. Los elegí de manera aleatoria, todos fueron guardados en archivos txt. Los  discursos varían, tanto en cantidad de palabras como el tipo de evento en el cual se emitió.

La muestra de discursos fue tratada tanto para hacer un Corpus Global, como uno por cada Mes, y de igual forma se compararon todos los discursos de manera independiente para identificar similaridad.

Lo que decidí hacer.

Escribí unas funciones las cuales me permitían dos cosas, extraer las palabras que aparecen con mayor frecuencia y aplicando LDA tomé 20 palabras del primer tópico detectado. Definí un corpus por año,  por mes y tomé los discursos qué más escuche comentar, el discurso del Informe de Gobierno y el Discursos en la ONU como la pareja inicial a comparar.

Así que use la medida de similaridad de Jaccard, que es una medida muy sencilla para definir entre los discursos cuales eran más cercanos según las palabras más frecuentes y cuales por las 20 palabras detectadas en el primer tópico.

Primero muestro un ejemplo y al final muestro los resultados obtenidos con todos los discursos.

Informe de Gobierno e informe ONU.

El código lo comparto más adelante. Los resultados son los siguientes:

#Comparación entre los discursos de la UNO y el Informe de Gobierno

dir1="C:\\.....ruta.....\\150928_ONU.txt"

dir2="D:\\.....ruta......\\150902_Informe_EPN.txt"

#Extracción del mensaje
Doc1<-msg(dir1)
Doc2<-msg(dir2)

#Contrucción de la matriz de terminos
TMD1<-tdm2(Doc1)
TMD2<-tdm2(Doc2)

#Se contruye su matrix de frecuencias
L1=TablaFreq(TMD1)
L2=TablaFreq(TMD2)

head(L1,5)
#    Términos  Frecuencia density ocurrencia
#356 naciones   15       0.01764706 1
#351 mundo      10       0.01176471 1
#330 méxico      9       0.01058824 1
#391 paz         9       0.01058824 1
#551 unidas      9       0.01058824 1
 
head(L2,5)
#     Términos Frecuencia density ocurrencia
#1597  méxico    92    0.014212884 1
#1603  mil       69    0.010659663 1
#1661  nacional  51    0.007878882 1
#1767  país      51    0.007878882 1
#1596  mexicanos 48    0.007415418 1

Jaccard(head(L1,20),head(L2,20))
#0.1428
#Gráfica
graphFreq(L1)
graphFreq(L2)

#Contruyo el DTM de cada texto
DTM_1=DTM(Doc1)
DTM_2=DTM(Doc2)

#Extraego las priemras 20 palabras asociadas con el primer tópico
top1<-TopicLDA(DTM_1)
top2<-TopicLDA(DTM_2)

top1
[1] "naciones" "derechos" "onu" "organización" "respeto" "agenda" 
[7] "frente" "humanidad" "humanos" "armas" "colectiva" "con" 
[13] "general" "seguridad" "señores" "acción" "clara" "consejo" 
[19] "favor" "futuro" 

top2
[1] "méxico" "con" "mil" "país" "reforma" 
[6] "gobierno" "mayor" "año" "administración" "familias" 
[11] "república" "años" "por" "justicia" "programa" 
[16] "educación" "crecimiento" "condiciones" "partir" "pobreza

Jaccard(top1,top2)
#0.0256

Se observa que hay una similaridad mayor entre los discursos cuando se consideran la palabras con mayor frecuencia o apariciones. La palabra que generó mayor crítica en los medios fue “populismo” y “populistas”. Estas dos palabras hicieron que se criticaran mucho los 2 discursos, lo que muestran los datos es que esas palabras aparecen muy poco, comparado con el efecto que generaron en los medios. De cierto modo son como “palabras singulares” que ejercieron mucho impacto en la apreciación del discurso.

La gráficas de las 50 palabras más citadas son las siguientes:

EPN_ONU

Discurso ONU

EPN_IG

Discurso Informe de Gobierno

Se observa que entre los dos discursos existe una diferencia considerable entre el tipo de palabras que se emplean y la frecuencia, también de la gráfica se puede apreciar que hay cierta “forma” diferente entre las palabras más citadas. Lo cual es notorio que la palabra “méxico” domina el discurso de Informe de Gobierno, por otro lado las palabras “naciones” y “mundo” el discurso de la ONU.

Lo que no muestra gran cambio es el desvanecimiento del color, que representa los cambios de la “densidad” de la frecuencia de las palabras. Concluir algo de estos datos es sutil y posiblemente confuso y atrevido, así que limitándome a lo estadístico se puede decir que hay en estructura diferencias pero muy sutiles, casi no se puede decir nada de esta comparación.

Dándole una interpretación a los tópicos y las mediciones de similaridad, los discursos en cuanto a “estructura”; las palabras más frecuentes, muestran un 14% de similaridad. Pero cuando uno analiza las palabras asociadas al primer tópico detectado, se observa que realmente solo son similares en escaso 2%. Los tópicos me parecen más relevantes, y es notoria la diferencia detectada. Las 20 palabras del primer tópico muestran al discurso de la ONU como algo “global” o “mundial”, y al discurso del informe de gobierno lo muestran como algo “nacional” y de problemas “socio económicos”.

En conclusión; como es de esperar,  se puede decir que los discursos no son tan similares (bajo esta medida de similaridad).

Lo global

Haciendo un corpus con los 48 discursos analizar, puedo comparar con respecto a los otros dos y analizar el comportamiento de la similaridad.

#Procesamiento de los discursos
#Librería para la nube de palabras
library(wordcloud)

dir="D:.....ruta....."
setwd(dir)
#Lista de directorios
filesall<-dir()
#Documentos y corpus
corpusalldoc<-sapply(filesall,function(p)msg(p))
corpusall<-tdm2(corpusalldoc)

#Tabla de frecuencias
Lall=TablaFreq(corpusall)

#Gráfica de frecuencias
graphFreq(Lall)

#Tópicos
DTMall<-DTM(corpusalldoc)
Topic_all<-TopicLDA(DTMall)

#Selección de conjunto de palabras para la nube de palabras
Grupo1<-rowSums(as.matrix(corpusall))
Grupo1<-subset(Grupo1,Grupo1>=30)
m<-as.matrix(Grupo1)
dim(m)
word_freqs = sort(rowSums(m), decreasing=TRUE)
dm = data.frame(word=names(word_freqs), freq=word_freqs)
wordcloud(dm$word, dm$freq, random.order=FALSE, random.color=FALSE,colors=brewer.pal(10,"Dark2"))

La gráfica y la nube de palabras que se obtiene es la siguiente:

Palabras_48 disc

Nube_48disc

Comparando con la métrica de similaridad se tienen lo siguiente:

#Medida de Similaridad
#Palabras más frecuentes
#UNO vs los 48 Discursos
Jaccard(Lall[,1],L1[,1])
#0.123
#Informe de Gobierno vs 48 Discursos
Jaccard(Lall[,1],L2[,1])
#0.282

#Tópicos
#ONU vs los 48 Discursos
Jaccard(Topic_all,top1)
#0.0526
#Informe de Gobiernos vs los 48 Discursos
Jaccard(Topic_all,top2)
#0.1428

Se aprecia que la medida de similaridad es mayor entre el discurso de informe que el de la ONU con respecto al corpus de los discursos. Era de esperar que fuera así, por la naturaleza del discurso y la fecha a la cual corresponde.

Es claro que los dos discursos tienen mayor similaridad con la muestra de discursos en el año, que entre ellos. Las dos preguntas que me hago al observar esto es, ¿cómo se comporta esta medida de similaridad por mes? y ¿cuál discurso muestra mayor similaridad con el de gobierno y el de la ONU?

Por Mes y por Discurso.

Haciendo la comparación por meses, se tienen que a comparar los dos discursos, el de la ONU y el del Informe, se tienen que por corpus construido por mes, se tiene gráficas como las siguientes:

Frec_Meses

El comportamiento por tópicos muestra otro comportamiento, el cual genera la siguiente gráfica:

Topic_Meses

Las gráficas de la métrica por mes muestra lo que uno puede esperar, el discurso de la ONU llega a no tener nada de similaridad en los meses de Marzo y Abril con la frecuencia de palabras, pero peor aún muestra poca similaridad con el primer grupo de tópicos en los meses de Octubre, Diciembre, Enero, Abril, Mayo y Septiembre. Es decir, el discurso dado en la ONU considerando que su primer grupo de tópicos se refiere aspecto mundiales o globales, ese no fue tema en esos meses con respecto a la muestra.

Por otro lado, el discurso del informe uno espera que sea similar al dado cada año o que las palabras que se usan en el mes de Septiembre suelen ser usuales. Eso muestra la primera gráfica, pero además vuelve a ser similar al inicio del año. Por otro lado al considerar los tópicos no muestra el mismo comportamiento, resulta que el mes de Junio es por alguna razón el mes con el cual muestra mayor similaridad. Eso me resulta raro, pero así resultó la medida de similaridad.

Ahora considerando cada uno de los discursos elegidos para analizar, el comportamiento que muestran con la frecuencia de palabras es el siguiente:

Freq_por_discurso

Esta muestra algo un poco más interesante, primero el mes de Septiembre muestra mayor similaridad y me resulta extraño que solo con el mes de Septiembre en el 2014, pero también uno puede observar que el mes de Enero con respecto al informe muestra un comportamiento de alta similaridad. Uno puede pensar que con motivo de inicio del año los discursos suelen ser “alentadores” , “nacionalistas”, “de mejoras” , etc. Esto pienso que puede ser interesante revisar una muestra de varios años y comparar como se comporta conforme pasan los años y quizás muestra estacionalidad la medida de similaridad.

Respecto al informe de la ONU muestra que no es usual que en los discursos se haga uso del mismo tipo de palabras, lo cual uno puede esperarlo ya que no suele decirse mucho del contexto “global”, como en el discursos de la ONU.

La gráfica de los tópicos, muestra el siguiente comportamiento:

 Topic_por_discursoLos tópicos muestran una cosa curiosa, el discurso de la ONU muestra entre los meses de Noviembre-Diciembre una alza, ¿efecto de la navidad para hablar del mundo?..no lo se, esto igual permite que si se hace una muestra mayor analizar si hay algún efecto en el primer grupo de tópicos detectados con la técnicas LDA.

Por otro lado el comportamiento del discurso del Informe muestra una alta similaridad en meses como Enero, Mayo y Junio. De nuevo que el mes de Enero aparezca con valores considerablemente mayores me hace suponer que al inicio del año y a medio año suelen tener este comportamiento de realzar o motivar ciertas cosas “nacionalistas” o “de esperanza” de ser un mejor país. No lo se, solo me hace pensar que teniendo una muestra mayor uno puede empezar hacer cosas más interesantes y jugar a poner algunas hipótesis para experimentar.

 ¿Qué cosas hacer para mejorar esto?

Haciéndome auto-críticas, pienso que hacer una muestra mayor y con discursos de varios años puede resultar más interesante. Por otro lado hacer uso de mejores técnicas o de otras técnicas de medidas de similaridad para explorar como se comportan los discursos con varias medidas. Por último no estaría mal hacer una muestra de otros mandatarios para revisar como evolucionan los tópicos y ver como se comportan ante camios o hechos históricos, cosas de ese estilo.

 Código

Comparto las funciones principales, el resto son muchas líneas de código de loops o de procesar un poco los datos para hacer las gráficas. Por lo cual comparto solo lo más importante del código.

library(tm)
library(NLP)
library(ggplot2)

#######################################
#Mensaje
msg<-function(path){
 con<-file(path,open='rt')
 text<-readLines(con, warn = FALSE, n=-1, encoding = "UCS-2LE")
 close(con)
 return(paste(text, collapse = "\n"))
}

#####################################
#TDM

tdm<-function(doc){
 control<-list(removeWords(stopwords("spanish")),
 removePunctuation,
 removeNumbers,
 content_transformer(tolower),
 minDocFreq=2)
 doc.corpus<-Corpus(VectorSource(doc))
 doc.tdm<-TermDocumentMatrix(doc.corpus,control)
 return(doc.tdm)
}

#######################################
#TDM versión 2
tdm2<-function(doc){
 docCor<-Corpus(VectorSource(doc))
 docs <- tm_map(docCor, stripWhitespace)
 docs <- tm_map(docs, removeWords, stopwords("spanish"))
 docs <- tm_map(docs, removePunctuation)
 docs <-tm_map(docs,removeNumbers)
 docs <- tm_map(docs,content_transformer(tolower))
 DocsTDM <- TermDocumentMatrix(docs) 
 return(DocsTDM)
 }

############################################
#Tabla de frecuencias
TablaFreq<-function(TDM){
 docmatrix <- as.matrix(TDM)
 doc.counts <- rowSums(docmatrix)
 doc.df <- data.frame(cbind(names(doc.counts),as.numeric(doc.counts)),stringsAsFactors = FALSE)
 names(doc.df) <- c("Términos", "Frecuencia")
 doc.df$Frecuencia <- as.numeric(doc.df$Frecuencia)
 doc.occurrence <- sapply(1:nrow(docmatrix),
 function(i)
 {
 length(which(docmatrix[i, ] > 0)) / ncol(docmatrix)
 })
 doc.density <- doc.df$Frecuencia / sum(doc.df$Frecuencia)
 doc.df <- transform(doc.df,density = doc.density,ocurrencia =doc.occurrence)
 S=head(doc.df[with(doc.df, order(-Frecuencia)),], n=50)
 return(S)
 }

##############################################
#Gráfica de frecuencias
graphFreq<-function(L){
 library(ggplot2)
 #Se debe de introducir la matriz con frecuencias
 #Gráfica de palabras y frencia
 p<-ggplot(L,aes(x=factor(Términos, levels=unique(as.character(Términos)) ), y=Frecuencia))+geom_bar(stat = "identity",aes(fill=density))+
 coord_flip()+xlab('Apariciones en el texto')+ylab('Las 50 palabras más frecuentes') 
 return(p)
 }

##########################################
#Función para extraer el Document term Matrix
DTM<-function(Texto){
 docCor<-Corpus(VectorSource(Texto))
 docs <- tm_map(docCor, stripWhitespace)
 docs <- tm_map(docs, removeWords, stopwords("spanish"))
 docs <- tm_map(docs, removePunctuation)
 docs<-tm_map(docs,removeNumbers)
 docs <- tm_map(docs,content_transformer(tolower))
 DocsTDM <- DocumentTermMatrix(docs) 
 return(DocsTDM)
}

################################################
#Topic
TopicLDA<-function(DTMdoc){
 #Regresa las 20 palabras relevantes del primer tópico detectado
 library(topicmodels)
 r.lda=LDA(DTMdoc,method="Gibbs",2)
 L=terms(r.lda,20)
 return(L[,1])
 }

##############################################
#Similaridad de Jaccard
Jaccard<-function(A,B){
 a=length(intersect(A,B))
 b=length(union(A,B))
 a/b
}

Análisis de Grafos

En esta entrada hago una recopilación de varios post de Andrie de Vries. Ya que el ejemplo que realiza se puede replicar sin problema alguno, la muestra de datos es suficientemente grande como para ser interesante y combina varias técnicas para detectar comunidades o clusters, con la finalidad de encontrar los nodos más sobre salientes del grafo.

Una área de investigación importante tanto en matemáticas, en física como en ingeniería; es la investigación en Grafos Grandes o redes largas (Large Graph). Sobre este tema recomiendo la referncia [1]. En lo personal es una área que me gusta y que pese a saber poco, trato de leer y de aprender al respecto, quizás más aspectos teóricos que prácticos; pero suele ser grato tener algunos ejemplos de donde aparecen los objetos matemáticos estudiados.

No puedo, ni pretendo hablar de teoría de grafos o gráficas por que no soy el adecuado para hablar a detalle del tema, pero recomiendo las referencias [2,3] para acercarse al tema desde el lado de algoritmos y desde el punto de vista de matemáticas discretas. Un libro en línea que puede ser una grata introducción a manejar grafos en python se puede ver en la referencia [4].

De que trata el ejemplo.

Hago uso de las idea y parte del código de Andrine Vries para mostrar como procesar un grafo. Este está formado por una muestra de paquetes o librerías de R, así que con más  7000 nodos es interesante aplicar ciertos algoritmos.

Las etapas son las siguientes:

  1. Cargar los datos y explorar parte de ellos.
  2. Aplicar el algoritmos PageRank.
  3. Aplicar técnicas de Cluster para grafos y comparar entre algunos algoritmos.
  4. Mostrar algunas representaciones gráficas sobre el grafo y el impacto de los algoritmos en él.

La idea es explicar el uso de cada librería en los puntos anteriores y el uso de Gephi.

Etapa 1. Generar los datos

Antes de mostrar un ejemplo sencillo de como generar un grafo, la idea de grafo es contar con nodos y con aristas. Donde las aristas conectan a los nodos o se quedan conectando al mismo nodo. Una idea gráfica se puede ver como la siguiente imagen:

usa-7-2

Grafo de libro Fractional Graph Theory

Pero también puede ser algo mucho más extraño, por la cantidad de nodos y la conexión entre las aristas. De modo tal que se pueden ver como:

map-of-internet

Visualización de la Intenet

Entonces la conexión entre redes y grafos pues es natural, prácticamente son lo mismo; pero la perspectiva con la cual la estudian desde el punto de vista de física es distinta a la que se estudia en matemáticas.

En R se tiene la librería miniCRAN con la cual se pueden extraer datos de los paquetes o librerías. Cuando se instala una librería muestra en la mayoría de los casos la solicitud de instalar las librerías con las cuales tiene dependencia en caso de no contar con ellas. Así que hay muchas librerías que tienen dependencia con muchas otras librerías y otras, simplemente con unas cuantas. La mayoría usan las librerías base de R, pero con miniCRAN se pueden ignorar estas y analizar solo las que no son básicas.

#Ejemplo 1
 librería(miniCRAN)

tags<-"ggplot2"

#Se extraen datos de las librerías con las cuales tiene relación de tres modos diferente, considerando o descartando algunas.

pkgDep(tags, suggests=FALSE, enhances=FALSE, includeBasePkgs = TRUE)

pkgDep(tags, suggests=FALSE, enhances=FALSE)

pkgDep(tags, suggests=TRUE, enhances=TRUE)

#Gráfica de las relaciones entre el paquete ggplot2

set.seed(1)

plot(makeDepGraph(tags, includeBasePkgs=FALSE, suggests=TRUE, enhances=TRUE), 
 legendPosEdge = c(-1, 1), legendPosVertex = c(1, 1), vertex.size=30, cex=0.7)

La gráfica que se obtiene es la siguiente:

ggplot2_graphs

Si se replica el código pero para una lista de librerías, se puede pensar en buscar el grafo de ciertas librerías que se usan para tareas relacionadas, ejemplo dplyr y reshape2, ggplot2 y lattice; etc. Pero para mostrar que no todos los paquetes están relacionados se puede considerar en la lista algunas librerías que tienen un uso diferente, como tm y lars.

#Ejemplo 2
library(miniCRAN)

tags<-c("dplyr","reshape2","ggplot2","tm","lattice","data.table","lars")

set.seed(1)

plot(makeDepGraph(tags, includeBasePkgs=FALSE, suggests=TRUE, enhances=TRUE), 
 legendPosEdge = c(-1, 1), legendPosVertex = c(1, 1), vertex.size=15,cex=0.5)

La gráfica que se obtienen es la siguiente:

dplyr-reshape2-data.table-ggplot2-lattice-tm_graphs

Observando la imagen se observa que hay dos nodos; tm y lars, que se muestran “disjuntos” a la red que muestran los otros paquetes. En este caso, se aprecia gráficamente que entre las librerías para procesar datos y generar gráficas existen ciertas relaciones, las cuales pueden ser analizadas.

Entonces la idea es extraer el grafo de todos los paquetes de R hasta hoy 29-09-2015 y aplicar algunos algoritmos para jugar con él. La cantidad de paquetes disponibles es de 7234 hasta hoy, lo cual hace pensar en que una imagen gráfica de dicho grafo es algo complicada para detectar cuales paquetes se aglomeran más que otros, o cuales son casi asilados del resto (como en las imágenes anteriores).

La extracción de los datos se realiza con el siguiente código:

library(miniCRAN)
library(magrittr)

MRAN <- "http://mran.revolutionanalytics.com/snapshot/2015-09-29/"

#Se extrae la matriz con todos los paquetes
pdb <- MRAN %>%
 contrib.url(type = "source") %>%
 available.packages(type="source", filters = NULL)

#Para visualizar el tipo de dato y las primeras 5 filas de la matriz de datos
class(pdb)
head(pdb)

#Se construye un grafo con la matriz

g <- pdb[, "Package"] %>%
 makeDepGraph(availPkgs = pdb, suggests=FALSE, enhances=TRUE, includeBasePkgs = FALSE)

#Se puede ver el tipo de dato
class(g)
head(g)

Etapa 2.Aplicación de PageRank

Por el tamaño del grafo resulta conveniente aplicar algún tipo de algoritmo para detectar algo sobre la relación entre paquetes. La idea de aplicar el algoritmo PageRank para detectar los paquetes que tienen más relevancia en la red.

Para ver detalles sobre el algoritmo les recomiendo la referencia [5]. En breve la idea del algoritmo es asignar un valor numérico a cada paquete y por medio de eso poder definir cuales tienen mayor relevancia según el valor que se le asignó que es un modo de sintetizar el peso o impacto del paquete en los otros paquetes. Esto no le hace el merito adecuado a dicho algoritmo, pero los detalles teóricos se salen de poder explicarse en la entrada.

Aplicando el algoritmo que se tienen implementado en la librería igraph, se obtienen los siguientes resultados:

#Aplicación de pageRank
library(igraph)

#Se aplica el algoritmo
pr <- g %>%
 page.rank(directed = FALSE) %>%
 use_series("vector") %>%
 sort(decreasing = TRUE) %>%
 as.matrix %>%
 set_colnames("page.rank")

head(pr, 10)
# Los 10 paquetes más importantes
#   page.rank
#Rcpp 0.021333171
#MASS 0.019112827
#Matrix 0.009390048
#ggplot2 0.009136545
#mvtnorm 0.008061911
#lattice 0.007902018
#survival 0.007814334
#plyr 0.007025619
#sp 0.004832534
#XML 0.004669243

Haciendo una representación gráfica del orden implementado por el algoritmo PageRank, se obtienen lo siguiente:

#Top 10
set.seed(40)
pr %>%
 head(10) %>%
 rownames %>%
 makeDepGraph(pdb) %>%
 plot(main="Top de paquetes", cex=0.5)

#Top 30
set.seed(40)
pr %>%
 head(30) %>%
 rownames %>%
 makeDepGraph(pdb) %>%
 plot(main="Top 30 de paquetes", cex=0.5)

Las gráficas que se obtienen son las siguientes:

top_10

Top_30

Se observa que el grafo de los primeros 10 es más o menos claro, pero que el de los primeros 30 ya empieza a mostrar cierta complejidad o se muestra más complicado para determinar visualmente algunos aspectos, pese a eso se aprecian en la imagen un par de paquetes asilados.

Etapa 3.-Detectando comunidades.

Los clusters, son parte de las técnicas de aprendizaje no supervidado en Machine Learning. Existen muchas técnicas o algoritmos, algunos al paso del tiempo ya sea por la eficiencia computacional o por el tipo de resultados que han permitido obtener terminan siendo predominantes o persistir ante la aparición de nuevas técnicas o algoritmos.

En la entrada “Análisis de Cluster, un ejemplo sencillo” mostré como usar dicha técnica con unos conjuntos de datos. La situación cambia un poco cuando se trabaja sobre grafos o redes, ya que se tiene en particular un tipo de datos (el grafo) sobre el cual se hace uso de sus propiedades para definir los algoritmos para detectar clusters.

En las referencias [6,7,8,9] se puede leer varios exploraciones experimentales sobre el tipo de algoritmos que más predominan en las investigaciones y al implementar dichos algoritmos. Aplico dos algoritmos para detectar comunidades o clusters en el grafo y mido la similaridad entre las dos comunidades con la medida de similaridad de Jaccard para comparar los resultados obtenidos entre los dos algoritmos.

Siguiendo a Andrie de Vries, solo tomo el 80% de la red de los cuales se eligen los nodos más importantes según el algoritmo PageRank.

#Se aplica al mismo grafo ordenado por PageRank pk
#Se selecciona el 80% de los datos

cutoff <- quantile(pr[, "page.rank"], probs = 0.2)
popular <- pr[pr[, "page.rank"] >= cutoff, ] 
toKeep <- names(popular)

vids <- V(g)[toKeep]
gs <- induced.subgraph(g, vids = toKeep)

#Aplicación del primer algoritmo walktrap community
#Se obtienen 2118 comunidades

cl <- walktrap.community(gs, steps = 3)

#Se ordenas
topClusters <- table(cl$membership) %>% 
 sort(decreasing = TRUE) %>% 
 head(25)

#Se gráfica el comportamiento de los cluster o de las comunidades

plot(topClusters, main="Tamaño Cluster", ylab="Numero de cluster", type="b", lwd=2)

#Se aplica el segundo algoritmo, infoMaps community
#Se obtienen 2280 comunidades

cl1 <- infomap.community(gs)

topClusters1 <- table(cl1$membership) %>% 
 sort(decreasing = TRUE) %>% 
 head(25)

plot(topClusters1, main="Tamaño de Clusters", ylab="Numero de clusters", type="b", lwd=2)

Cluster_walktrap

InfoMaps_clusters

Las gráficas muestran el comportamiento de los cluster o comunidades, lo cual muestra que se pueden elegir los primero 10 clusters como los más importantes. Siguiendo a Andrie de Vries se eligen los 10 clusters más importantes.

Para comparar el comportamiento entre los dos algoritmos; además de que el primero detecta 2118 comunidades y el segundo 2280, comparo los 10 cluster más importantes por medio del índice de similaridad de Jaccard. Para los que tienen nociones de teoría de conjuntos, la idea es muy sencilla es comparar cuantos elementos están en la intersección dividido entre los elementos en la unión. Esto da una medida de similaridad, con la cual comparo entre las 10 comunidades más importantes por algoritmo.

#Eligiendo los 10 clusters más importantes

cluster <- function(i, clusters, pagerank, n=10){
 group <- clusters$names[clusters$membership == i]
 pagerank[group, ] %>% sort(decreasing = TRUE) %>% head(n)
}

#Para walktrap
z1 <- lapply(names(topClusters)[1:10], cluster, clusters=cl, pagerank=pr, n=20)

#Para infoMaps
z2 <- lapply(names(topClusters1)[1:10], cluster, clusters=cl1, pagerank=pr, n=20)

#Comparación entre los cluster de cada algoritmo

Jaccard<-function(A,B){
 a=length(intersect(A,B))
 b=length(union(A,B))
 a/b
 }
s=vector(length=10)

for(i in 1:10){
 L1=names(z1[[i]])
 L2=names(z2[[i]])
 s[i]=Jaccard(L1,L2)
 }

library(ggplot2)
t=1:10
qplot(x=t,y=s,geom=c("line", "point"), main="Similaridad entre los dos algoritmos", xlab="Número de Cluster",ylab="Valor de Similaridad")

Similaridad

Se observa que en los dos algoritmos en el primer cluster se tienen más similaridad, la lista de librerías por algoritmo son las siguientes:

#Para el algoritmo Walktrap
names(z1[[1]])
 [1] "MASS" "Matrix" "mvtnorm" "lattice" 
 [5] "survival" "igraph" "coda" "nlme" 
 [9] "rgl" "boot" "RColorBrewer" "numDeriv" 
[13] "Hmisc" "mgcv" "cluster" "gtools" 
[17] "car" "fields" "lme4" "xtable" 

#Para el algoritmo infoMaps
names(z2[[1]])
 [1] "Rcpp" "plyr" "XML" "stringr" 
 [5] "reshape2" "foreach" "jsonlite" "rJava" 
 [9] "cluster" "car" "fields" "lme4" 
[13] "corpcor" "e1071" "doParallel" "vegan" 
[17] "quadprog" "rpart" "randomForest" "plotrix"

La interpretación o explicación de los resultados puede varias, ejemplo en los conjuntos de librerías anteriores aparece una combinación entre librerías para procesar datos, para graficar, para aplicar modelos predictivos o de clasificación, etc.

Para visualizar el grafo de las librerías, la intención era mostrarlo con Gephi pero tuve un problema con el programa y hasta que después repare o vea que es lo que tiene comparto la visualización de toda la red. Pero con igraph se puede hacer una imagen buena de como se comportan los cluster en la red.

#Uso de igraph, parte de código de Andrei de Vries

#Se carga la librería
library(igraph)

# Se preparan los datos para construir la gráfica

cl$degree <-(degree(gs)[cl$names])
cl$cluster<-unname(ave(cl$degree, cl$membership,FUN=function(x)names(x)[which.max(x)]))
V(gs)$name <- cl$cluster

E(gs)$weight <- 1
V(gs)$weight <- 1

gcon<-contract.vertices(gs, cl$membership,vertex.attr.comb = list(weight = "sum", name = function(x)x[1], "ignore"))

gcon<-simplify(gcon, edge.attr.comb = list(weight = "sum", function(x)length(x)))

gcc<-induced.subgraph(gcon, V(gcon)$weight > 20)

V(gcc)$degree<-unname(degree(gcc))

#construcción de la gráfica
par(mar = rep(0.1, 4)) 
g.layout<-layout.kamada.kawai(gcc)
plot.igraph(gcc, edge.arrow.size = 0.1, layout = g.layout, vertex.size = 0.5 * (V(gcc)$degree))

La imagen de grafo que se obtiene es el siguiente:

Grafo1_vis

Se muestra en la imagen como la librería con círculos mayores tienen mayor relevancia en la parte de la red, por su puesto que no están las 7000 librerías representadas pero esta imagen da idea de cuales son las librerías más importantes o las más usadas, partiendo de la clasificación hecha con el algoritmo walktrap y pagerank.

Espero la entrada de una idea como tratar en general datos que pueden ser analizados como un grafo, cabe mencionar que cuando se analiza un corpus de textos o un texto, este también puede ser analizado y explorado visualizando su matriz de términos como la matriz asociada a un grafo y eso hace que se tengan más herramientas para analizar lo que sucede con dicho texto o la relevancia de sus tópicos.

Referencias:

  1. Large Network and graph limits.Lásló Lovasz
  2. Estructura de datos y algoritmos.
  3. Graph Theory
  4. http://interactivepython.org/runestone/static/pythonds/index.html
  5. http://www.mmds.org/
  6. http://www-personal.umich.edu/~mejn/papers/epjb.pdf
  7. http://arxiv.org/abs/1004.3539
  8. http://deim.urv.cat/~clara.granell/publications/ijcss_clara.pdf
  9. http://www.syssec-project.eu/m/page-media/3/moradi-sea12.pdf

Regresión en Python

La regresión.

En varias entradas de Series de Tiempo y de Machine Learning en R comenté sobre la regresión lineal. El método es bien conocido en varios campos y quizás es el más natural de aprender, ya que la idea de trazar la recta que mejor describe la relación entre puntos de un plano parece familiar e intuitivo.

Ejemplo_opt_regresion

El problema con el método o el modo en que se enseña pude variar, como el método en ocasiones se subvalora debido a que ese considera en general que es solo “agregar línea de tendencia” a los datos, pero se deja de lado la revisión de los estadísticos y de revisar el comportamiento de los errores de la estimación.

En la entrada de manera breve comento un pequeño recorrido sobre las técnicas de regresión clásica y las penalizadas. Las principales referencias sobre las cuales comento los métodos son [1,2] y algunos artículos o reportes que dejo en las referencias. Lo que ya no comento es la parte estadística; es decir, los detalles sobre las pruebas de hipótesis y las potencias de las pruebas las comenté en otras entradas. En esta trato de solo mencionar lo mínimo sobre esos detalles, pero enfocarme más en otros usos y las variaciones a la técnica tradicional.

El modelo en general.

La idea de la regresión es tener un conjunto de datos de entrada y explicar o describir las salidas como una combinación lineal de los datos de entrada.

Es decir, se tienen los datos: Datos_entrada

y el modelo que se pretende describir es:

Modelo_Lineal

En la ecuación anterior se interpreta como la salida del modelo f(X) sea lo más parecida a Y, que son los valores reales a estimar.

Para el caso de los mínimos cuadrados se tiene que la función a minimizar es:

LQ

LQ2

Lo que se hace es sustituir la expresión que de f(x) que depende de las β’s para minimizar la ecuación, lo que se hace es tratar de encontrar los mejores β’s para describir Y.

En los métodos modernos lo que se penalizar la estimación de los β’s, debido a que al estimar la regresión por mínimos cuadrados sobre un subconjunto de entrenamiento ( train set) de los datos, es decir algunos de los x’s, permite obtener los parámetros  β’s sobre estimados  y al momento de usar el modelo para el conjunto de datos de prueba (test set) se obtienen tasas de errores mayores.

Los métodos con penalización usuales son Ridge y Lasso, los dos agregan un parámetro λ que implica el problema de tener un método para definir el mejor valor para el modelo. Explico con ejemplos que se hace para determinar el parámetro.

Para el modelo Ridge a la función de errores es como la de los mínimos cuadrados pero agregando un término más. La ecuación es la siguiente:

RR

Para el modelo Lasso es similar a la de Ridge, pero en la penalización en lugar de considerar los β’s  al cuadrado se considera su valor absoluto. La ecuación es la siguiente:

LR

Lamento poner las ecuaciones, pero creo que pese a que no se entienda la teoría de manera total , al observando las diferencias en las ecuaciones se puede apreciar la diferencia en los métodos. En la primera gráfica de la entrada no lo mencioné, pero si se presiona sobre la imagen y se observa con cuidado que tiene trazadas 3 rectas sobre los datos de color negro, estas tres rectas corresponden a los 3 métodos distintos.

Una visión o idea gráfica de lo que pasa con el modelo o técnica, es visto en la siguiente imagen. Al final el objetivo de la regresión lineal es encontrar el plano que mejor satisface ciertas condiciones para asegurar que predice o explica lo mejor posible los datos.

Regression_Lineal

Imagen tomada del libro de Trevor Hastie y Robert Tibshirani.

Ejemplo 1

Para visualizar la diferencia en las estimaciones haciendo uso de Python, sobre todo de la librería scikit-learn hago uso del código de Jaques Grobler, modifico parte de su código para agregar la estimación de la regresión Rigde y Lasso, además elijo una muestra mayor de datos para prueba para que la gráfica sea notoriamente más atractiva.

Lo que hace el código:

  1. Se toma una base de datos sobre diabetes la cual se puede cargar directamente en la librería scikit-learn.
  2. Se divide la muestra en datos de entrenamiento( train) y de prueba (test).
  3. Se estiman los 3 modelos de regresión y se estiman los estadísticos base de las estimaciones.
  4. Se hace la gráfica de los datos de test.
 # -*- coding: 850 -*-
print(__doc__)


# Code source: Jaques Grobler
# License: BSD 3 clause

#Se cargan los modulos

import matplotlib.pyplot as plt
import numpy as np
from sklearn import datasets, linear_model

# Se cargan los datos
diabetes = datasets.load_diabetes()


# Se usa cierta muestra de los datos
diabetes_X = diabetes.data[:, np.newaxis]
diabetes_X_temp = diabetes_X[:, :, 2]

# Se separa la muestra de datos
diabetes_X_train = diabetes_X_temp[:-30]
diabetes_X_test = diabetes_X_temp[-30:]

# Se separa las variables explicadas
diabetes_y_train = diabetes.target[:-30]
diabetes_y_test = diabetes.target[-30:]

# Se crean los modelos lineales
regr = linear_model.LinearRegression()
regr2=linear_model.Lasso(alpha=.5)
regr3=linear_model.Ridge(alpha=.5)

#Se entrenan los modelos
regr.fit(diabetes_X_train, diabetes_y_train)
regr2.fit(diabetes_X_train, diabetes_y_train)
regr3.fit(diabetes_X_train, diabetes_y_train)


print u'Regresión Mínimos Cuadrados Ordinarios'
# Coeficiente
print'Coeficientes:',regr.coef_
# MSE
print("Residual sum of squares: %.2f"
 % np.mean((regr.predict(diabetes_X_test) - diabetes_y_test) ** 2))
# Varianza explicada
print('Varianza explicada: %.2f\n' % regr.score(diabetes_X_test, diabetes_y_test))

print u'Regresión Lasso' 
# Coeficiente
print'Coeficientes:', regr2.coef_
# MSE
print("Residual sum of squares: %.2f"
 % np.mean((regr2.predict(diabetes_X_test) - diabetes_y_test) ** 2))
# Varianza explicada
print('Varianza explicada: %.2f\n' % regr2.score(diabetes_X_test, diabetes_y_test))

print u'Regresión Ridge'
#Coeficiente
print'Coeficientes:', regr3.coef_
# MSE
print("Residual sum of squares: %.2f"
 % np.mean((regr3.predict(diabetes_X_test) - diabetes_y_test) ** 2))
# Varianza explicada
print('Varianza explicada: %.2f\n' % regr3.score(diabetes_X_test, diabetes_y_test))


# Plot 
plt.scatter(diabetes_X_test, diabetes_y_test, color='black')
plt.plot(diabetes_X_test, regr.predict(diabetes_X_test), color='blue',
 linewidth=3, label=u'Regresión MCO')
plt.plot(diabetes_X_test, regr2.predict(diabetes_X_test), color='yellow',
 linewidth=3, label=u'Regresión Lasso')
plt.plot(diabetes_X_test, regr3.predict(diabetes_X_test), color='green',
 linewidth=3, label=u'Regresión Ridge')
plt.title(u'Regresión Lineal por 3 metodos diferentes')
plt.legend()
plt.xticks(())
plt.yticks(())

plt.show()

Los resultados serían los siguientes:

Regresión Mínimos Cuadrados Ordinarios
Coeficientes: [ 941.43097333]
Residual sum of squares: 3035.06
Varianza explicada: 0.41

Regresión Lasso
Coeficientes: [ 720.32135439]
Residual sum of squares: 3265.36
Varianza explicada: 0.37

Regresión Ridge
Coeficientes: [ 612.64202818]
Residual sum of squares: 3458.17
Varianza explicada: 0.33

La gráfica se vería de la siguiente forma:

Tres_Métodos

El ejemplo anterior viendo la gráfica y los estadísticos tenemos como conclusión que le método MCO; mínimo cuadrados ordinarios, es el que mejor se ajusta a nuestros datos. En este caso uno se pregunta sobre la ventaja de usar otro método para la regresión lineal, si los resultados son favorables para la técnica más conocida.

Lo que no se menciona en el código, o lo que se queda sin explicar es la elección de λ para los modelos Ridge y Lasso. Lo que se hace para aplicar un buen modelo de estos dos métodos, es estrenar el modelo de modo tal que el valor de λ sea el “mejor”. Para esto se usa Cross-Validation, así que haciendo ligeros cambios al código anterior se puede mejorar la aplicación de los modelos lineales.

 # -*- coding: 850 -*-
print(__doc__)


# Code source: Jaques Grobler
# License: BSD 3 clause

#Se cargan los módulos
import matplotlib.pyplot as plt
import numpy as np
from sklearn import datasets, linear_model

# Se cargan los datos
diabetes = datasets.load_diabetes()


# Se eligen las muestras de datos
diabetes_X = diabetes.data[:, np.newaxis]
diabetes_X_temp = diabetes_X[:, :, 2]

# Se dividen los datos en training/testing 
diabetes_X_train = diabetes_X_temp[:-30]
diabetes_X_test = diabetes_X_temp[-30:]

# Se divide la variable explicada en  training/testing 
diabetes_y_train = diabetes.target[:-30]
diabetes_y_test = diabetes.target[-30:]

# Se crean los modelos lineales
regr = linear_model.LinearRegression()
#regr2=linear_model.LassoCV(alpha=)
regr3=linear_model.RidgeCV(alphas=[0.1,0.2,0.5,1.0,3.0,5.0,10.0])

# Se entrenan los modelos
regr.fit(diabetes_X_train, diabetes_y_train)
regr2=linear_model.LassoCV(cv=20).fit(diabetes_X_train, diabetes_y_train)
regr3.fit(diabetes_X_train, diabetes_y_train)


print u'Regresión Mínimos Cuadrados Ordinarios'
# Coeficiente
print'Coeficientes:',regr.coef_
# MSE
print("Residual sum of squares: %.2f"
 % np.mean((regr.predict(diabetes_X_test) - diabetes_y_test) ** 2))
# Varianza Explicada
print('Varianza explicada: %.2f\n' % regr.score(diabetes_X_test, diabetes_y_test))

print u'Regresión Lasso' 
# Coeficientes
print'Coeficientes:', regr2.coef_
# MSE
print("Residual sum of squares: %.2f"
 % np.mean((regr2.predict(diabetes_X_test) - diabetes_y_test) ** 2))
# Varianza Explicada
print('Varianza explicada: %.2f\n' % regr2.score(diabetes_X_test, diabetes_y_test))

print u'Regresión Ridge'
# Coeficiente
print'Coeficientes:', regr3.coef_
# MSE
print("Residual sum of squares: %.2f"
 % np.mean((regr3.predict(diabetes_X_test) - diabetes_y_test) ** 2))
# Varianza Explicada
print('Varianza explicada: %.2f\n' % regr3.score(diabetes_X_test, diabetes_y_test))


# Plot 
plt.scatter(diabetes_X_test, diabetes_y_test, color='black')
plt.plot(diabetes_X_test, regr.predict(diabetes_X_test), color='blue',
 linewidth=3, label=u'Regresión MCO')
plt.plot(diabetes_X_test, regr2.predict(diabetes_X_test), color='yellow',
 linewidth=3, label=u'Regresión Lasso')
plt.plot(diabetes_X_test, regr3.predict(diabetes_X_test), color='green',
 linewidth=3, label=u'Regresión Ridge')
plt.title(u'Regresión Lineal por 3 metodos diferentes')
plt.legend()
plt.xticks(())
plt.yticks(())

plt.show()

Lo que se tienen como resultado es lo siguiente:

Regresión Mínimos Cuadrados Ordinarios
Coeficientes: [ 941.43097333]
Residual sum of squares: 3035.06
Varianza explicada: 0.41

Regresión Lasso
Coeficientes: [ 940.48954236]
Residual sum of squares: 3035.57
Varianza explicada: 0.41

Regresión Ridge
Coeficientes: [ 850.17738251]
Residual sum of squares: 3103.11
Varianza explicada: 0.40

La gráfica correspondiente es:

Tres_Métodos2

 

Lo que se aprecia es que los modelos mejoran considerablemente, si bien el modelo MCO sigue siendo uno de los mejores para este ejemplo, algo a considerar es que es una muestra pequeña y de pocas dimensiones. Así que lo interesante es aplicar la regresión a muestras mayores y de datos con una mayor cantidad de variables.

Ejemplo 2-Regresión Multiple

Para el siguiente ejemplo hago uso de una matriz de datos de tamaño 506 líneas de datos con 12 variables de entrada y una variable que se trata de estimar por métodos lineales. Los datos corresponden a los valores de las viviendas en los suburbios de Boston.

Tomo estos datos para mostrar como aplicar la regresión a datos con varia variables, para este caso no hago una gráfica pero se pueden revisar los datos o valores de los estadísticos para tener un modo de elegir entre modelos.

El código es muy similar al anterior y lo que se puede hacer como variación es aplicar alguna técnica de reducción de dimensiones para elegir las variables que “impactan realmente” en el modelo o jugar con ejecutar la regresión no con todas las variables, sino empezar con unas y poco a poco agregar nuevas para revisar como se comporta el modelo al estimar con solo unas cuantas variables. Algo que tampoco hago en este caso es explorar las variables antes de aplicar la regresión, pero la intención es ver que prácticamente es igual la ejecución de una regresión multi lineal.

#import matplotlib.pyplot as plt
import numpy as np
from sklearn import datasets, linear_model

#13 variables explicativas

# Se cargan los datos
Boston = datasets.load_boston()



# Se divide en training/testing sets
Boston_train = Boston.data[:405,]
Boston_test = Boston.data[405:,]


# La variable explicativa se divide en training/testing sets
Boston_y_train = Boston.target[:405]
Boston_y_test = Boston.target[405:]

regr = linear_model.LinearRegression()
regr2=linear_model.LassoCV(cv=20).fit(Boston_train, Boston_y_train)
regr3=linear_model.RidgeCV(alphas=[0.1,0.2,0.5,0.7,1.0,1.5,3.0,5.0,7.0,10.0])
#np.size(Boston_y_train = Boston.target[:405])
regr.fit(Boston_train, Boston_y_train)
regr3.fit(Boston_train,Boston_y_train)

print u'Regresión Mínimos Cuadrados Ordinarios'
#Coeficiente
print'Coeficientes:',regr.coef_
# MSE 
print("Residual sum of squares: %.2f"
 % np.mean((regr.predict(Boston_test) - Boston_y_test) ** 2))
# Varianza explicada
print('Varianza explicada: %.2f\n' % regr.score(Boston_train, Boston_y_train))


print u'Regresión Ridge'
# Coeficientes
print'Coeficientes:', regr3.coef_
# MSE
print("Residual sum of squares: %.2f"
 % np.mean((regr3.predict(Boston_test) - Boston_y_test) ** 2))
# Varianza Explicada
print('Varianza explicada: %.2f\n' % regr3.score(Boston_train, Boston_y_train))


print u'Regresión Lasso' 
# Coeficiente
print'Coeficientes:', regr2.coef_
# MSE
print("Residual sum of squares: %.2f"
 % np.mean((regr2.predict(Boston_test) - Boston_y_test) ** 2))
# Varianza Explicada
print('Varianza explicada: %.2f\n' % regr2.score(Boston_train, Boston_y_train))

Los resultados del código son los siguientes resultados:

Regresión Mínimos Cuadrados Ordinarios
Coeficientes: [ -1.94651664e-01 4.40677436e-02 5.21447706e-02 1.88823450e+00
 -1.49475195e+01 4.76119492e+00 2.62339333e-03 -1.30091291e+00
 4.60230476e-01 -1.55731325e-02 -8.11248033e-01 -2.18154708e-03
 -5.31513940e-01]
Residual sum of squares: 33.58
Varianza explicada: 0.74

Regresión Ridge
Coeficientes: [ -1.93961816e-01 4.42902027e-02 4.66864588e-02 1.87777854e+00
 -1.37158180e+01 4.76760163e+00 1.65803755e-03 -1.28477100e+00
 4.57255634e-01 -1.57415769e-02 -7.97849246e-01 -1.78511420e-03
 -5.33113498e-01]
Residual sum of squares: 33.15
Varianza explicada: 0.74

Regresión Lasso
Coeficientes: [-0.13537099 0.04758408 -0. 0. -0. 3.26851037
 0.01014632 -0.9131925 0.39587576 -0.01796467 -0.65864709 0.00511243
 -0.67828222]
Residual sum of squares: 23.48
Varianza explicada: 0.71

Si se tienen cuidado se puede observar que los dos primeros modelos muestran mucha “cercanía”, los dos tienen estimaciones casi similares y estadísticos muy parecidos. El modelo Lasso al contrario, parece considerablemente diferente y con valores considerablemente diferentes. El punto clave es cómo elegir entre los modelos. En esta parte sería adecuado hacer uso de algunos indicadores sobre la calidad del modelo como AIC o BIC, etc.

Ejemplo 3. Usando la regresión de otro modo.

El siguiente ejemplo es un ejemplo que considero bastante bonito y muestra el uso de Lasso  y Ridge que en general no se piensa, respecto el código y el reconocimiento a Emmanuelle Gouillart, sobre la cual recomiendo revisar su tutorial sobre procesamiento de imágenes. No hago cambio alguno al código, lo que hago es explicar que hace en cada etapa para clarificar y valorar las regresiones con penalización como herramientas de reconstrucción de imágenes.

El código:

print(__doc__)

# Author: Emmanuelle Gouillart <emmanuelle.gouillart@nsup.org>
# License: BSD 3 clause

import numpy as np
from scipy import sparse
from scipy import ndimage
from sklearn.linear_model import Lasso
from sklearn.linear_model import Ridge
import matplotlib.pyplot as plt


def _weights(x, dx=1, orig=0):
 x = np.ravel(x)
 floor_x = np.floor((x - orig) / dx)
 alpha = (x - orig - floor_x * dx) / dx
 return np.hstack((floor_x, floor_x + 1)), np.hstack((1 - alpha, alpha))


def _generate_center_coordinates(l_x):
 l_x = float(l_x)
 X, Y = np.mgrid[:l_x, :l_x]
 center = l_x / 2.
 X += 0.5 - center
 Y += 0.5 - center
 return X, Y


def build_projection_operator(l_x, n_dir):
 """ Compute the tomography design matrix.

 Parameters
 ----------

 l_x : int
 linear size of image array

 n_dir : int
 number of angles at which projections are acquired.

 Returns
 -------
 p : sparse matrix of shape (n_dir l_x, l_x**2)
 """
 X, Y = _generate_center_coordinates(l_x)
 angles = np.linspace(0, np.pi, n_dir, endpoint=False)
 data_inds, weights, camera_inds = [], [], []
 data_unravel_indices = np.arange(l_x ** 2)
 data_unravel_indices = np.hstack((data_unravel_indices,
 data_unravel_indices))
 for i, angle in enumerate(angles):
 Xrot = np.cos(angle) * X - np.sin(angle) * Y
 inds, w = _weights(Xrot, dx=1, orig=X.min())
 mask = np.logical_and(inds >= 0, inds < l_x)
 weights += list(w[mask])
 camera_inds += list(inds[mask] + i * l_x)
 data_inds += list(data_unravel_indices[mask])
 proj_operator = sparse.coo_matrix((weights, (camera_inds, data_inds)))
 return proj_operator


def generate_synthetic_data():
 """ Synthetic binary data """
 rs = np.random.RandomState(0)
 n_pts = 36.
 x, y = np.ogrid[0:l, 0:l]
 mask_outer = (x - l / 2) ** 2 + (y - l / 2) ** 2 < (l / 2) ** 2
 mask = np.zeros((l, l))
 points = l * rs.rand(2, n_pts)
 mask[(points[0]).astype(np.int), (points[1]).astype(np.int)] = 1
 mask = ndimage.gaussian_filter(mask, sigma=l / n_pts)
 res = np.logical_and(mask > mask.mean(), mask_outer)
 return res - ndimage.binary_erosion(res)


# Generate synthetic images, and projections
l = 128
proj_operator = build_projection_operator(l, l / 7.)
data = generate_synthetic_data()
proj = proj_operator * data.ravel()[:, np.newaxis]
proj += 0.15 * np.random.randn(*proj.shape)

# Reconstruction with L2 (Ridge) penalization
rgr_ridge = Ridge(alpha=0.2)
rgr_ridge.fit(proj_operator, proj.ravel())
rec_l2 = rgr_ridge.coef_.reshape(l, l)

# Reconstruction with L1 (Lasso) penalization
# the best value of alpha was determined using cross validation
# with LassoCV
rgr_lasso = Lasso(alpha=0.001)
rgr_lasso.fit(proj_operator, proj.ravel())
rec_l1 = rgr_lasso.coef_.reshape(l, l)

plt.figure(figsize=(8, 3.3))
plt.subplot(131)
plt.imshow(data, cmap=plt.cm.gray, interpolation='nearest')
plt.axis('off')
plt.title('original image')
plt.subplot(132)
plt.imshow(rec_l2, cmap=plt.cm.gray, interpolation='nearest')
plt.title('L2 penalization')
plt.axis('off')
plt.subplot(133)
plt.imshow(rec_l1, cmap=plt.cm.gray, interpolation='nearest')
plt.title('L1 penalization')
plt.axis('off')

plt.subplots_adjust(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0,
 right=1)

plt.show()

La imagen que se obtiene es la siguiente:

Img1_Regression

Primero explicando el resultado lo que se tienen es una imagen que se crea y se toma como imagen original, la imagen del centro y la derecha son aproximaciones lineales a la original. Lo interesante es que las dos aproximaciones son buenas, pero la generada por Lasso es mucho mejor que la generada por Ridge.

¿Cómo hace esto el código? En general la idea de aplicar una regresión es para predecir datos y se piensa en la imagen de la recta o el plano que se aproxima a los datos que se tratan de modelar.

En esta caso también se hace una aproximación lineal, lo cual es generado por las funciones:

1.-generate_synthetic_data()
2.-build_projection_operator()

Entonces en resumen, lo primero que se hace es cargan las librerías requeridas, lo cual es normal en todos los códigos. Luego se definen funciones que hacen uso de las funciones base de numpy para procesar y construir matrices. Explotando esto, lo que se construyen las dos funciones mencionadas anteriormente.

Así la muestra de datos queda considerada dentro de la variable proj sobre la cual se aplica el método np.ravel(), con lo cual se aplican las dos regresiones penalizadas para  hacer la aproximación lineal con la imagen original que se gurda en la variable data.

Lamento no explicar todo a mucho detalle, pero creo que los comentario de Emmanuelle son claros y los pocos conceptos que están ocultos en el código son claros con un poco de conocimiento en álgebra lineal.

Espero que aprecien el uso de los dos métodos de regresión lineal penalizados siendo usados en una aproximación lineal a un conjunto de datos diferente a lo que suele uno ver en los libros de texto.

Ejemplo 4.-Otro ejemplo con ciertos datos.

En el libro de Willi Richert y Luis Pedro Coelho, dedican dos capítulos a mostrar el uso de la regresión lineal como una herramienta para  predecir el  rating  en un sistema de recomendación. Recomiendo leer los capítulos 7 y 8 para que sea claro como se hace uso de la regresión para éste tipo de sistemas, referencia [3]. Como ejemplo de como usar la regresión es bastante bueno,  pero algo largo y no creo bueno compartirlo en el blog ya que son varios script para mostrar como funcionaria el sistema. Por ello recomiendo comprar el texto y revisar el código que comparten en su repositorio en Github para replicar y estudiar dicho ejemplo.

Un punto clave para los sistemas de recomendación es considerar el caso cuando se tienen una cantidad pequeña de usuario comparada con una gran cantidad de productos. Ejemplo, suponiendo que se va el supermercado, se tienen una cantidad de quizás 5 mil artículo o productos y se reciben en dicha tienda solo mil clientes. Construyendo la matriz con mil filas y 5 mil columnas, se tienen un problema del tipo Ρ»Ν.

Este tipo de problemas suelen aparecer mucho cuando uno analiza datos de los empleados de una empresa, sobre productos o cierto tipo de inventarios o matrices de datos que tienen que ver con los procesos de algún empresa.

En caso de considerar la regresión como un buen modelo aplicar a este tipo de datos, el problema es que se tiene una matriz de datos que no permitiría hacer una división de entre train/ test para el algoritmo, entonces como alternativa lo que se hace es emplear una técnica de Cross-Validation la cual permite elegir los parámetros de modelo con “mejor calidad”.

La técnica es llamada k-fold cross-validation, lo que se hace es que se divide la muestra en k particiones y se elige una de esas como el conjunto de test y las k-1 como train. Se repite el proceso k veces y al final se eligen las k-medias de los coeficientes como los valores de los coeficientes del modelo. Dejo en las referencias los tutoriales de Andrew Moore, en el cual pueden ver en las presentaciones como funciona de manera general el método.

CV-kfold

Ejemplo gráfico de k-fold Croos-Validation

Sobre el código, entonces como ejemplo hago uso del código de Willi Richert y Luis Pedro Coelho con una ligera modificación para considerar el uso de Lasso y comparar los resultados.

#Módulos
import numpy as np
from sklearn.metrics import mean_squared_error, r2_score
from sklearn.datasets import load_svmlight_file
from sklearn.linear_model import LinearRegression, LassoCV
from sklearn.cross_validation import KFold

#Cargar los datos
data, target = load_svmlight_file('data/E2006.train')

lr = LinearRegression()

# Scores:

#Estimación del modelo lineal
lr.fit(data, target)

#Estimación del modelo Lasso
lr2=LassoCV(cv=5).fit(data,target)

#Estimaciones o predicciones
pred = lr.predict(data)
pred2 = lr2.predict(data)


print('RMSE sobre los datos de entrenamiento, {:.2}'.format(np.sqrt(mean_squared_error(target, pred))))

print('R2 sobre los datos de entrenamiento, {:.2}'.format(r2_score(target, pred)))
print('')

pred = np.zeros_like(target)
kf = KFold(len(target), n_folds=5)
for train, test in kf:
 lr.fit(data[train], target[train])
 pred[test] = lr.predict(data[test])

print('RMSE sobre los datos de prueba (5-fold), {:.2}'.format(np.sqrt(mean_squared_error(target, pred))))
print('R2 sobre los datos de prueba (5-fold), {:.2}'.format(r2_score(target, pred)))

#Estimaciones para Lasso
print('RMSE en Lasso(5-fold), {:.2}'.format(np.sqrt(mean_squared_error(target, pred2))))
print('R2 en Lasso(5-fold), {:.2}'.format(r2_score(target, pred2)))

El código regresa cierta información, al cuales son los estadísticos básicos para saber como se comporta el modelo.

RMSE sobre los datos de entrenamiento, 0.0012
R2 sobre los datos de entrenamiento, 1.0

RMSE sobre los datos de prueba(5-fold), 0.59
R2 sobre los datos de prueba(5-fold), 0.12
RMSE en Lasso(5-fold), 0.37
R2 en Lasso(5-fold), 0.66

Si se revisan los primeros datos causa confusión, por que parecen que el primer modelo es “perfecto”, pero eso es una mala señal sobre el modelo. Para más detalles dejo algunas referencias sobre sobre los estadísticos y su valoración.

Entonces con este último ejemplo la intención es mostrar que cuando se tienen condiciones del tipo Ρ»Ν  lo adecuado es usar Cross-Validation para este tipo de muestras al momento de aplicar la regresión.

Varios errores que cometemos al aplicar una regresión lineal.

En un post de Cheng-Tau Chu listó 7 errores que pasan al aplicar algoritmos de Machine Learning. Eligiendo los errores que se presentan en la regresión lineal en general, la lista sería algo así:

  • No revisar la correlación de los errores
  • No tener varianza constante  en los errores.
  • No considerar la importancia de los outliers.
  • No cuidar los problemas de la multicolinealidad o colinealidad en los datos.
  • Considerar un modelo lineal cuando las interacciones de los datos son no-lineales.

Hablar de cada punto requiere dar un ejemplo de donde falla el modelo lineal y como abordar dicho problema. Recomiendo la referencia [2] para revisar algunos detalles.

Conclusión: Espero que los ejemplos den una idea general sobre como se usa la regresión lineal. Hay otras variaciones de la regresión que son LAR y Elastic Net, sobre las dos hay mucho material y sobre todo que son técnicas más nuevas y han estado en recientes cambios con nuevas investigaciones.

Referencias:

  1. https://mitpress.mit.edu/books/machine-learning-2
  2. http://statweb.stanford.edu/~tibs/ElemStatLearn/
  3. http://www.amazon.es/Building-Machine-Learning-Systems-Python/dp/1782161406
  4. http://www.springer.com/us/book/9780387310732

 

Suppor Vector Machine ( Máquina de Soporte de Vectores) en Python

Sobre SVM

Las técnicas de Máquina de Soporte de Vectores (SVM) por alguna extraña razón es uno de los algoritmos que más interés despiertan y sobre todo el que creo que genera más sobre estimación. Es posible que no he usado suficiente la familia de algoritmos como para reconocer tanto su valor.

Hace no más de un año un amigo en una propuesta sugería usar esta familia de algoritmos sin revisar si realmente convenía hacer uso de ellos para el tipo de datos que se analizarían. La justificación de este amigo era que SVM eran los “mejores” algoritmos, lo cual es falso, no se puede calificar a una familia de algoritmos como los mejores, se requiere probar varios contra los datos para determinar cuales funcionan mejor con ellos.

Dos ejemplos que he compartido donde involucro SVM se pueden revisar en las entradas:

1.-Comparacion-entre-svm-naive-bayes-arboles-de-decision-y-metodos-lineales.

2.-Maquina-de-soporte-vectorial-svm-sopport-vector-machine

En las dos entradas uso SVM para clasificar, si bien sus resultados son buenos en las muestras de datos no resulta ser el mejor método al compararlo con otra familia de algoritmos. Lo que recomiendo es probar con más de una familia de algoritmos para analizar alguno de los problemas y elegir más de un estadístico o indicador que puedan servir para comparar la eficiencia de los algoritmos.

Si se tiene interés en conocer más detalles sobre los algoritmos recomiendo consultar la referencia [1,2,3,4], en orden con respecto a la dificultad del contenido y de las explicaciones. Los primeros 3 ejemplos se pueden encontrar en la documentación de la librería scikit-learn.

La relación entre las técnicas o algoritmos es la siguiente; el primero es Margen Maximal de Clasificación (MMC), su generalización es Soporte de Vectores para Clasificación (SVC) y esta a su vez se generaliza a SVM.

El primer ejemplo con datos simulados, son la imagen clásica para explicar como funciona SVM, pero lo que se muestra en sí es el caso MMC.

El código es el siguiente:

#MMC 
print(__doc__)

#Librerías requeridas

import numpy as np
import matplotlib.pyplot as plt
from sklearn import svm

# Se crean los datos
np.random.seed(0)
X = np.r_[np.random.randn(70, 2) - [2, 2], np.random.randn(70, 2) + [2, 2]]
Y = [0] * 70 + [1] * 70

#Se estima el modelo
clf = svm.SVC(kernel='linear')
clf.fit(X, Y)

# Se construye la recta que separa las clases
w = clf.coef_[0]
a = -w[0] / w[1]
xx = np.linspace(-5, 5)
yy = a * xx - (clf.intercept_[0]) / w[1]


# support vectors
b = clf.support_vectors_[0]
yy_down = a * xx + (b[1] - a * b[0])
b = clf.support_vectors_[-1]
yy_up = a * xx + (b[1] - a * b[0])

#Gráfica
plt.plot(xx, yy, 'k-')
plt.plot(xx, yy_down, 'k--')
plt.plot(xx, yy_up, 'k--')

plt.scatter(clf.support_vectors_[:, 0], clf.support_vectors_[:, 1],
 s=80, facecolors='none')
plt.scatter(X[:, 0], X[:, 1], c=Y, cmap=plt.cm.Paired)

plt.title('Recta separadora')
plt.axis('tight')
plt.show()

La gráfica es la siguiente:

MMC

A lo que se refiere el algoritmo con el nombre de Margen Maximal, es al margen que se obtiene entre las dos clases de datos, siendo el margen la distancia entre las líneas punteadas.

Esto en teoría es una optimización sobre el margen y la generalización se da cuando no se puede obtener una línea recta que no intercepte o que separe de manera clara las muestras de datos. El ejemplo clásico se hace con los datos Iris, los cuales permiten mostrar como se generaliza el uso de los margenes que separan los datos, cuando no son líneas rectas.

El código muestra el ejemplo de como se separan de manera no lineal las fronteras.

#Código

print(__doc__)

#Módulos
import numpy as np
import matplotlib.pyplot as plt
from sklearn import svm, datasets

#Se importan los datos
iris = datasets.load_iris()
X = iris.data[:, :2] 
y = iris.target

h = .02

#Parámetro de regularización
C = 1.0 

svc = svm.SVC(kernel='linear', C=C).fit(X, y)
rbf_svc = svm.SVC(kernel='rbf', gamma=0.7, C=C).fit(X, y)
poly_svc = svm.SVC(kernel='poly', degree=3, C=C).fit(X, y)
lin_svc = svm.LinearSVC(C=C).fit(X, y)

#Probé pero no me gustó el resultado con sigmoid
#sigmoid_svc=svm.SVC(kernel='sigmoid').fit(X,y)

#Se crean los marcos para las gráficas

x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1
y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1
xx, yy = np.meshgrid(np.arange(x_min, x_max, h),np.arange(y_min, y_max, h))

#Títulos para las gráficas
titles = ['SVC con kernel lineal',
 'Lineal SVC',
 'SVC con RBF kernel',
 'SVC con polinomio(grado 3) kernel']


for i, clf in enumerate((svc,lin_svc, rbf_svc, poly_svc)):
 # Se grafican las fronteras 
 plt.subplot(2, 2, i + 1)
 plt.subplots_adjust(wspace=0.4, hspace=0.4)

 Z = clf.predict(np.c_[xx.ravel(), yy.ravel()])

 #Color en las gráficas
 Z = Z.reshape(xx.shape)
 plt.contourf(xx, yy, Z, cmap=plt.cm.Paired, alpha=0.8)

 #Puntos de entrenamiento
 plt.scatter(X[:, 0], X[:, 1], c=y, cmap=plt.cm.Paired)
 plt.xlabel('Longitud Sepal')
 plt.ylabel('Peso Sepal')
 plt.xlim(xx.min(), xx.max())
 plt.ylim(yy.min(), yy.max())
 plt.xticks(())
 plt.yticks(())
 plt.title(titles[i])

plt.show()

La gráfica que se obtiene es la siguiente:

SVC_Iris

El otro uso es para estimar la regresión y un ejemplo visual que se aprecia rápido es en la regresión de datos en dos dimensiones, es decir; en el plano.

El código es el siguiente:

#Generación de datos

print(__doc__)
import numpy as np
from sklearn.svm import SVR
import matplotlib.pyplot as plt

#Generación de datos
X = np.sort(5 * np.random.rand(100, 1), axis=0)
y = np.sin(X).ravel()
y[::5] += 3 * (0.5 - np.random.rand(20))

#Modelos regresión fit

svr_rbf = SVR(kernel='rbf', C=1e3, gamma=0.1)
svr_lin = SVR(kernel='linear', C=1e3)
svr_poly = SVR(kernel='poly', C=1e3, degree=2)
y_rbf = svr_rbf.fit(X, y).predict(X)
y_lin = svr_lin.fit(X, y).predict(X)
y_poly = svr_poly.fit(X, y).predict(X)

#Mirando los resultados

plt.scatter(X, y, c='k', label='datos')
plt.hold('on')
plt.plot(X, y_rbf, c='g', label='Modelo RBF ')
plt.plot(X, y_lin, c='r', label='Modelo Lineal')
plt.plot(X, y_poly, c='b', label='Modelo Polinomial')
plt.xlabel('datos')
plt.ylabel('target')
plt.title('M. De Soport de Vectores')
plt.legend()
plt.show()

La gráfica que se muestra es la siguiente:

SVM_Regresión

En las primeras ocasiones que uno usa esta familia de algoritmo la confusión se genera con el concepto de “kernel”,  de modo no formal y burdo la idea de lo que hace el kernel es modificar el modo en el cual se logran definir los planos separadores o las fronteras entre clases de datos. Esto en palabras puede no parecer claro, pero el ejemplo natural para ejemplificar esto es generar una muestra de datos altamente no lineales y aplicar SVM con diferentes kernel lo cual mostrará como se comporta el algoritmo al cambiar el tipo de kernel. En el ejemplo hago una comparación con el kernel rbf, lineal, polinomial de grado 4 y 5, y el kernel sigmoid.

El código es el mismo para las cuatro estimaciones solo se cambiar una línea de código la cual para obtener las gráficas se comentan las no requeridas. El código es prácticamente una ligera modificación al ejemplo que se puede encontrar en scikit-learn.

El código para realizar este ejemplo es el siguiente:

#Código
#Module
print(__doc__)

import numpy as np
import matplotlib.pyplot as plt
from sklearn.svm import SVC
import time

start_time=time.time()
xx, yy = np.meshgrid(np.linspace(-3, 3, 1000),
 np.linspace(-3, 3, 1000))
np.random.seed(0)
X = np.random.randn(1000, 2)
Y = np.logical_xor(X[:, 0] > 0, X[:, 1] > 0)


# fit the model

#clf =SVC(kernel='poly',degree=5)
clf =SVC(kernel='poly',degree=4)
#clf =SVC(kernel='linear')
#clf=SVC()
clf.fit(X, Y)

# plot the decision function for each datapoint on the grid
Z = clf.decision_function(np.c_[xx.ravel(), yy.ravel()])
Z = Z.reshape(xx.shape)

plt.imshow(Z, interpolation='nearest',
 extent=(xx.min(), xx.max(), yy.min(), yy.max()), aspect='auto',
 origin='lower', cmap=plt.cm.PuOr_r)
contours = plt.contour(xx, yy, Z, levels=[0], linewidths=2,
 linetypes='--')
plt.scatter(X[:, 0], X[:, 1], s=30, c=Y, cmap=plt.cm.Paired)
plt.xticks(())
plt.yticks(())
plt.axis([-3, 3, -3, 3])
plt.title('Kernel Polinomio grado 4')
plt.show()

print(" The time required is %f seconds "%(time.time()-start_time))

Las gráficas que se obtienen son las siguientes:

SVM_Kernel_Lineal SVM_Kernel_Polinomial_degree4
SVM_Kernel_RBF

SVM_Kernel_Polinomial_degree5

Lo que se observa en las 4 gráficas es que el modelo con kernel rbf y un kernel polinomial de grado 4, parecen ser los dos mejores modelos para clasificar los datos en dos clases.Se observa que el modelo lineal no es para nada apropiado y que le cuesta trabajo al SVM con ese kernel poder separar los datos. El único modo que considero o creo que sirve para detectar si los datos son altamente lineales o altamente no lineales es hacer una exploración en los datos.

Un par de ejemplos.

Tomando datos de la plataforma http://datos.gob.mx/,  hago un ejemplo burdo,  ya que los datos que tomo son sobre los delitos por estado o Entidad Federativa en distintas clasificaciones de delitos, estos datos cubren desde el mes de Enero del 2001 hasta diciembre del 2006.  Lo que hago es medir la media de todos los delitos por cada uno de los registros de los datos, después teniendo la media defino una escala simple para asignar una variable indicadora a los estados que tienen mayor y menos tasa delictiva.

Con lo anterior defino un dos clases, debido a que no uso los datos para hacer un  análisis cuidadoso sino solo para comparar el uso del algoritmo, entonces omito la relevancia de tener las fechas de la medición por mes y año, tomo una muestra aleatoria de aproximadamente el 70% para entrenar el algoritmo y el resto para probar la eficiencia.

En resumen las etapas son las siguientes:

  1. Procesar y limpiar los datos.
  2. Asignar la variable de clasificación o indicadora.
  3. Entrenar el algoritmo y comparar con el conjunto de prueba.
  4. Comparar contra otro algoritmo, Naive Bayes.

 Primero recomiendo descargar los datos, ya que se cuenta con ellos mediar la media por cada fila y teniendo esa tabla la cargamos a python.

#Procesamiento
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import os
from pandas import Series

#Cargamos los datos
os.chdir('Directorio donde se pusieron los datos')

Datos=pd.read_csv('Archivos_con_los_datos.csv',sep=',',decimal=',')

Datos.head()
Datos.shape

#Exploración

Datos['Media'].hist()
plt.xlabel('Valores de la Media de la Delincuencia por Estado')
plt.ylabel('Valor de la Media')
plt.title('Histograma')
plt.show()

#Limpiamos los datos por los Missing values

sum(Datos.isnull().values.ravel())
#31

#Cleaning of data 
Datos1=Datos.dropna()
Datos1.shape 
#(2273,27)

S=Datos1['Media']
#Exploración 2
S2=[log(y) for y in S]

plt.hist(S2,col='green')
plt.xlabel('Valores de la Media de la Delincuencia por Estado')
plt.ylabel('Valor de la Media')
plt.title('Histograma')
plt.show()

#Missing Values
sum(Datos.isnull().values.ravel())
#31

#Cleaning of data 
Datos1=Datos.dropna()
Datos1.shape 
#(2273,27)

Series(S).mean()
#9.32


Hist_Delic_Mex

Histograma de las Medias de la Delincuencia por Estado

Hist_Del_Mex_Log

Histograma de los Logaritmos de las Medias de la Delincuencia por Estado

Después de explorar de modo breve como se comportan los datos ahora medimos a media creamos una nueva variable y le asignamos un valor dependiendo de la medición de la media. Agrego también los valores del logaritmo de la media ya que esta variable muestra que el comportamiento de las medias tienen una distribución de colas pesadas.

#Nueva variable
Datos1['Log_Media']=S2

Datos1.head()
Datos1.var()
 
#total data
 
len(Datos1[Datos1['Media']>9.3169]['Log_Media'])
#569
Datos1['Indicador']=0
Datos1.ix[Datos1.Media>9.3169,'Indicador']=1

Datos1.shape
#(2273,29)

Ahora preparo los datos para entrenar el algoritmo y aplico SVM, tanto con un kernel lineal como con uno no lineal.

# Preparación de datos y SVM
from sklearn.svm import SVC
from sklearn import svm
from sklearn.naive_bayes import GaussianNB
import random 

rows=random.sample(Datos1.index,1536)

Datos_train=Datos1.ix[rows]
Datos_test=Datos1.drop(rows)

type(Datos_train)
X=np.asarray(Datos_train.icol(range(4,26)))
Y=np.asarray(Datos_train.icol(28))
X_test=np.asarray(Datos_test.icol(range(4,26)))
Y_test=np.asarray(Datos_test.icol(28))
X_test

#Aplico los algoritmos de SVM

#Kernel no lineal
clf=SVC()
clf.fit(X,Y)
len(Y_test)
#737
P=clf.predict(X_test)
sum(P==Y_test)/737.0
#0.778

#Kernel Lineal
lin_clf=svm.LinearSVC()
lin_clf.fit(X,Y)
Pred=lin_clf.predict(X_test)
sum(Pred==Y_test)/737.0
#0.9701

#Naive Bayes
gnb = GaussianNB()
Prediction=gnb.fit(X,Y).predict(X_test)
(Prediction!=Y_test).sum()
(Prediction==Y_test).sum()
683.0/737.0
#0.926

Se observa que el porcentaje de eficiencia de los algoritmos es el siguiente, para SVM con un kernel no lineal es del 77.8%, para SVM con un kernel lineal es del 97.01% y para el algoritmo Naive Bayes es de 92.6%.

Entonces en resumen, para esta muestra de datos se obtienen que SVM con un kernel lineal es el más eficiente como método de clasificación.

La moraleja de siempre es usar más de un algoritmo para una muestra de datos, ya que ninguno es el “mejor” solamente es más adecuada para cada cierto tipo de datos.

Referencias:

1.- An Introduction to Statistical Learning

2.-The Elements of Statistical Learning

3.-Machine Learning: A Probabilistic Perspective

4.-Pattern Recognition  and Machine Learning

Cursos y Referencias de Deep Learning

Decidí compartir información de lo que he estado revisando de Deep Learning, que va desde cursos, libros, manuales, ejemplos.

Pienso que uno debe de acercarse a estos temas ya que si bien hasta el momento están implementadas en sistemas grandes, como empresas financieras, Google, Facebook, etc. Pensando en lo que ha pasado con Machine Learning, que ahora resulta fácil hacer un ejemplo con una máquina estándar, lo que  supongo es que pronto ese tipo de situaciones serán viables para hacer uso de técnicas e infraestructuras para proyectos pequeños o medianos de Deep Learning.

Entonces tratando de poner un granito de arena en este mar de conocimiento me propuse compartir los cursos que he estado revisando y las referencias que están disponibles en internet.

En vídeo

Hay muchas platicas, desde introductorias hasta avanzadas sobre Deep Learning, para mi gusto el curso completo, con buenas referencias y explicaciones es el de Nando de Freitas de la universidad de Oxford.

El curso completo se puede seguir en sus canal de youtube:

https://www.youtube.com/user/ProfNandoDF

Una segunda opción son los cursos de Hugo Larochelle, lo único que recomiendo es tener paciencia ya que son vídeos de menor duración pero suelen contener más detalles técnicas.

Su curso de puede ver e su canal de youtube

https://www.youtube.com/user/hugolarochelle/featured

y sobre sus trabajos y cursos se puede ver en su página:

http://www.dmi.usherb.ca/~larocheh/index_en.html

Notas respecto al tema

Respecto a notas, textos y conocimientos necesarias para abordar el tema, lo que suelo hacer es revisar lo que indica alguna de las universidades top o la líder en el tema, en esta caso recomiendo seguir los consejos de la Universidad de Montreal que es la líder en este campo.

La ventaja es que todo está en la web, lo cual permite ser consultado a cualquier hora y de manera libre. En algunos casos se requiere un poco de formación en Matemáticas para entender algunos detalles de los algoritmos y en otros un poco de familiaridad con respecto a Python, para manejar Theano, pero las recomendaciones son las siguientes:

Libro

El mejor libro que hasta ahora he encontrado que hable del tema está disponible para ser leído en red, pero no para ser descargado en algún formato es el de Y. Bengio.

Existen capítulos en algunos muy buenos libros como el de Kevin P. Murphy, el cual es creo un breve resumen para tener una idea general antes de abordar la lectura de un texto más amplio.

Mi recomendación es el capítulo 28 del libro de Kevin P. Murphy

http://mitpress-ebooks.mit.edu/product/machine-learning

Y el mejor libro de Deep Learning , para mi gusto, lo pueden leer en la siguiente liga:

http://www.iro.umontreal.ca/~bengioy/dlbook/

Cursos y otras referencias

Los pocos cursos que he visto disponibles son los de la Universidad de Stanford, entonces aprovechando que los dejaron para revisar el material de manera libre sobre decir que se deben de aprovechar.

Los más interesantes a mi gusto son los siguientes dos:

Deep Learing para Procesamiento de Lenguaje Natural

Convolución de Redes Neuronales para reconocimiento visual

Uno menos recientes es:

Deep Learing para Procesamiento de Lenguaje Natural( sin magia)

Algunos blog o referencias que son bastante buenas son las siguientes:

Notas de Michael Nilsen

Tutorial Aprendizaje no Supervizado y Deep Learning

El mejor lugar para revisar un amplio material sobre Deep Learing, es decir; desde aspectos prácticos hasta nivel investigación están en el siguiente sitio( es una belleza):

http://deeplearning.net/reading-list/

Comentario final

Espero que con el material anterior se tenga una buena referencia sobre el tema, mi intención es poner en esta categoría algunos ejemplo de como usar los algoritmos y aclarar de modo breve cual es la idea central. Trataré de hacer una entrada por tema del curso de Nando de Freitas usando Theano y PyLearn2.

Mi única intensión es invitar acercarse al tema y practicar, ya que datos para hacerlo sobra, pero lo falta en ocasiones es el tiempo y la dedicación ha estar aprendiendo sobre el tema.

Comparación entre Máquina de Soporte Vectorial, Naive Bayes, Árboles de Decisión y Métodos Lineales

De lo que trata esta entrada.

En la entrada no explico a detalles técnicas sobre cada una de los algoritmos empleados, pero pueden consultarse en las referencias o en las categorías Machine Learning en R project y Machine Learning en Python.

Para el ejemplo solo uso código en R, la intención es usar tres muestras de datos simulados para mostrar las estimaciones que realiza la función svm de la librería e1071, solo en el caso de usarla para clasificar. Y por último uso unos datos  de correos clasificados como spam o no-spam provenientes del repositorio UCI Machine Learning Repository los cuales se encuentran cargados en la librería kernlab de  R project para aplicar las técnicas y comparar svm con otras técnicas.

Ejemplo.-Solo Suppor Vector Machine

En este ejemplo solo hago dos muestras de datos simulados, con distribuciones  Gaussianas y la intención es compartir los detalles de qué es en rasgos generales lo que hace el algoritmo de maquina de soporte vectorial al clasificar datos, esto es para dos dimensiones. Es decir, sobre un plano ya que visualmente se vuelve claro y con esa idea se puede pensar en el tipo de cosas que hace el algoritmo en más dimensiones.

Los primeros datos no se mezclan, por lo cual es visualmente claro como separarlos.

#Datos
#Librerías requeridas para el ejemplo

library(ggplot2)
library(e1071)

#Datos de mezcla 1

N<-500
x1<- rnorm(N) * 0.1
y1<- rnorm(N) * 0.1
X1<- t(as.matrix(rbind(x1, y1)))
x2<- rnorm(N) * 0.1 + 0.5
y2<- rnorm(N) * 0.1 + 0.5
X2<- t(as.matrix(rbind(x2, y2)))
X<- as.matrix(rbind(X1, X2))
X=data.frame(X)
X['Clus']=1
X[500:1000,3]=2
ggplot(X,aes(x=x1,y=y1))+geom_point(size=3,aes(colour=factor(Clus)))+ggtitle('Mezcla 1')+theme(plot.title=element_text(lineheight = 2,face='bold'))+
 xlab('Variable X1')+ylab('Variable Y1')

La gráfica de los puntos se ve así:Mezcla1_Gas

Para la siguiente mezcla considero tres muestras donde las clasifico en dos categorías.

#Mezcla 2
N<-500
x1<- rnorm(500,mean=0,sd=3)*0.1 
y1<- rnorm(500,mean=0,sd=3)*0.1 
X1<- t(as.matrix(rbind(x1, y1)))
x3<- rnorm(250,mean=0,sd=3)*0.1 
y3<- rnorm(250,mean=0,sd=3)*0.1+0.7 
X3<- t(as.matrix(rbind(x3, y3)))
x2<- rnorm(N,mean=1.5,sd=4)*0.1 + 0.5
y2<- rnorm(N,mean=0,sd=3)*0.1 + 0.5
X2<- t(as.matrix(rbind(x2, y2)))
Datos2<- as.matrix(rbind(X1, X3,X2))
Datos2=data.frame(Datos2)
Datos2['Clus']=1
Datos2[600:1250,3]=2
ggplot(Datos2,aes(x=x1,y=y1,colour=factor(Clus)))+geom_point(size=3)+ggtitle('Mezcla 2')+theme(plot.title=element_text(lineheight = 2,face='bold'))+
 xlab('Variable X1')+ylab('Variable Y1')

Mez2_GaussLa gráfica de los datos se ve como la anterior imagen.

Los datos muestran comportamientos distintos, los primeros es claro que se puede separar por una recta, pero la gráfica de la segunda muestra no se ve claramente que lo mejor sea separar los datos por una recta.

Como una buena práctica se debe de tomar una muestra de datos para entrenar el algoritmo (train set) y una cantidad de datos de prueba(test set). Para los siguientes ejemplos considero todos los datos, para ilustrar cual sería la gráfica de los datos predichos por el modelo.

La técnica de SVM para clasificar tienen de fondo la idea de encontrar el mejor “hiperplano” por medio del cual separar los datos. En este ejemplo el concepto de hiperplano es una recta, es decir; los datos están en dos dimensiones (x,y) y se busca la mejor recta que separe los datos, si uno piensa en datos con tres dimensiones (x,y,z) que pueden pensarse como alto, largo y ancho lo que se busca es encontrar el mejor plano que separa los datos, por ello el nombre de hiperplano.

Uso la librería e1071 para estimar el algoritmo SVM con el kernel lineal y radial para la primera muestra de datos, para la segunda uso lineal, polinomial, radial y sigmoid. La intención de este ultimo es compara cual de las cuatro muestra una similitud gráfica más parecida a los datos originales y comparo la tasa de predicciones correctas. Existen más librerías para hacer uso de SVM, dejo en la referencia las ligas.

Primera muestra de datos.

#SVM con Kernel Lineal
#Se procede a dejar una semilla para hacer una elección del mejor valor de parámetro cost

set.seed (1)
#Se hace cross-validation de k-fold, con k=10. Esto mediante la función tune()

tune.out=tune(svm,factor(Clus)~.,data=Datos1,kernel="linear",type='C-classification',scale=FALSE,ranges=list(cost=c(0.001,0.01, 0.1, 1,5,10,100)))
tune.out
summary(tune.out)

#Se elige el mejor modelo con el mejor valor para el parámetro cost

bestmod=tune.out$best.model
summary(bestmod)
# Gráfica del mejor modelo con kernel lineal

plot(bestmod,Datos1)

#Predicción
Pred=predict(bestmod,Datos1)

table(predicción=Pred,Valores_reales=Datos1$Clus)

# Valores_reales
# predicción 1 2
#         1 499 1
#         2 0 500
#Con los valores obtenidos se tiene un 99.9% de eficiencia 

Algunas observaciones, el modelo pasa por un proceso de validación cruzada usando la función tune(), esto siempre es recomendable para elegir un buen modelo. El parámetro cost, significa el nivel de penalización que permite el modelo, en otras palabras el algoritmo busca la mejor recta que separe los datos y como tal el costo para encontrar esa recta requiere tolerar posibles datos que afectan a la estimación de dicha recta. Esto quizás es una explicación burda y mala, pero detrás de todo algoritmo de Machine Learning está un proceso de optimización el cual determina el valor de los parámetros que hacen que el algoritmo tenga el mejor valor posible.

La gráfica que se obtiene de este modelo es la siguiente:

SVM_Muestra1_ker-lineal

 

Esta gráfica muestra las dos clases que se buscaban definir.La eficiencia es muy alta, de 1000 datos clasifica correctamente el 99.9%.

Lo que se espera de otro kernel es que prevalezca la eficiencia de SVM y más aún que la gráfica de la clasificación sea casi igual a la obtenido cuando se usa un kernel lineal.

#Kernel radial
set.seed (1)
tune.out=tune(svm,factor(Clus)~.,data=Datos1,kernel="radial",type='C-classification',scale=FALSE,ranges=list(cost=c(0.001,0.01, 0.1, 1,5,10,100)))
tune.out
summary(tune.out)
bestmod_r=tune.out$best.model
summary(bestmod_r)
plot(bestmod_r,Datos1)
Pred=predict(bestmod_r,Datos1)
table(predicción=Pred,Valores_reales=Datos1$Clus)
# Valores_reales
# predicción 1 2
# 1 499 1
# 2 0 500
#Eficiencia de la clasificación 99.9%

Al obtener la gráfica del modelo por SVM se tiene:

SVM_Muestra1_ker-radial

Lo cual muestra una imagen muy similar a la obtenida con el kernel lineal. La tasa de eficiencia al clasificar es prácticamente la misma, 99.9%

Con la segunda muestra de datos pensar en separar por una recta no parece lo natural. Primero hago la estimación con cuatro tipo de kernels de la función svm y muestro la gráfica que regresa el modelo.

#Segunda muestra de datos
#Se usa primero el kernel lineal

set.seed (1)
tune.out=tune(svm,factor(Clus)~.,data=Datos2,kernel="linear",type='C-classification',scale=FALSE,ranges=list(cost=c(0.001,0.01, 0.1, 1,5,10,100)))
tune.out
summary(tune.out)
bestmod_l=tune.out$best.model
summary(bestmod_l)
plot(bestmod_l,Datos2)
Pred=predict(bestmod_l,Datos2)
table(predicción=Pred,Valores_reales=Datos2$Clus)
# Valores_reales
# predicción 1 2
# 1 504 99
# 2 95 552
#Eficiencia de la clasificación 84.4%

Se tiene la siguiente gráfica como resultado de implementar el algoritmo:

SVM_Mezcla2_ker-lineal

Si se compara los datos que predice el modelo con respecto a los originales se tienen lo siguiente:

Orig_vs_Ker-lineal

Se observa que la predicción muestra totalmente separadas las dos clases y cabe notar que se tiene el 84.4% eficiencia.

#Kernel polinomial

set.seed (1)
tune.out=tune(svm,factor(Clus)~.,data=Datos2,kernel="polynomial",type='C-classification',scale=FALSE,ranges=list(cost=c(0.001,0.01, 0.1, 1,5,10,100)))
tune.out
summary(tune.out)
bestmod_p=tune.out$best.model
summary(bestmod_p)
plot(bestmod_p,Datos2)
Pred=predict(bestmod_p,Datos2)
table(predicción=Pred,Valores_reales=Datos2$Clus)
# Valores_reales
# predicción 1 2
# 0 543 174
# 1 56 477
#Eficiencia de la clasificación 81.6%, con kernel polinomial

La gráfica que se obtiene del modelo es:

SVM_Mezcla2_ker-polin

Y otra vez haciendo una gráfica para comparar la predicción con los datos originales se tiene:

Orig_vs_Ker-pol

Se aprecia que no es tan definida la recta que separa las clases como en el ejemplo del kernel lineal, pero la eficiencia se reduce ya que es de 81.6%

#Kernel sigmoid

set.seed (1)
tune.out=tune(svm,factor(Clus)~.,data=Datos2,kernel="sigmoid",type='C-classification',scale=FALSE,ranges=list(cost=c(0.001,0.01, 0.1, 1,5,10,100)))
tune.out
summary(tune.out)
bestmod_s=tune.out$best.model
summary(bestmod_s)
#Informa del tipo de modelo, del que puede considerarse como el mejor modelo
plot(bestmod_s,Datos2)
Pred=predict(bestmod_s,Datos2)
table(predicción=Pred,Valores_reales=Datos2$Clus)
# Valores_reales
# predicción 1 2
#         1 509 102
#         2 90 549
#Eficiencia de la clasificación 84.6%, con kernel Sigmoid

La gráfica que se obtiene del modelo es la siguiente:

SVM_Mezcla2_ker-sigmoid

Haciendo la gráfica comparativa de predicción y datos originales se obtiene:Orig_vs_Ker-sigmoid

 

Se aprecia que es muy similar a la que se obtiene con el kernel lineal y hasta el nivel de eficiencia resulta muy aproximado, ya que es del 84.6%

#Kernel Radial
set.seed (1)
tune.out=tune(svm,factor(Clus)~.,data=Datos2,kernel="radial",type='C-classification',scale=FALSE,ranges=list(cost=c(0.001,0.01, 0.1, 1,5,10,100)))
tune.out
summary(tune.out)
bestmod_r=tune.out$best.model
summary(bestmod_r)
plot(bestmod_r,Datos2)
Pred=predict(bestmod_r,Datos2)
table(predicción=Pred,Valores_reales=Datos2$Clus)
# Valores_reales
# predicción 1 2
# 1 497 94
# 2 102 557
#Eficiencia de la clasificación 84.3%, con kernel lineal

La gráfica obtenida con este kernel es la siguiente:

SVM_Mezcla2_ker-radial

 

Y al comparar las predicciones con los datos originales se tiene:

Orig_vs_Ker-radial

Se aprecia que es similar a la obtenida por el kernel lineal y sigmoid, más aún la eficiencia resulta ser del 84.3%.

Entonces en resumen los datos al ser clasificados mediante SVM con distintos kernel resultó en estos datos resulta tener mejor eficiencia el kernel sigmoid, por un porcentaje mínimo sobre el lineal.

Para mostrar como se comporta SVM con una muestra de datos altamente no lineales o que resulta difícil separar por medio de una recta , genero una muestra más.

#Datos altamente no lineales

#Distribución uniforme 1250 valores
x1=runif(1250)-0.5
x2=runif (1250)-0.5
#Variable Indicadora
y=1*(x1^2-x2^2> 0)

#Gráfica donde se muestran los puntos pintados por etiqueta de y
#Se contruye un data.frame con los datos
X=data.frame(x1,x2,y)
ggplot(data=X,aes(x=x1,y=x2))+geom_point(aes(colour=factor(y),shape=factor(y)))+ggtitle('Datos Originales')+theme(plot.title=element_text(lineheight = 2,face='bold'))+
 xlab('Variable X1')+ylab('Variable X2')

La gráfica de los datos es la siguiente:

Datos_No-lineales

En esta muestra de datos se aprecia que separar las clases por una recta no parece lo más eficiente, para constatarlo estimo la clasificación con el kernel lineal.

#Kernel lineal
set.seed (1)
tune.out=tune(svm,factor(y)~.,data=X,kernel="linear",type='C-classification',scale=FALSE,ranges=list(cost=c(0.001,0.01, 0.1, 1,5,10,100)))
tune.out
summary(tune.out)
bestmod=tune.out$best.model
summary(bestmod)
#Informa del tipo de modelo, del que puede considerarse como el mejor modelo
plot(bestmod,X)
Pred=predict(bestmod,X)
table(predicción=Pred,Valores_reales=X$y)
# Valores_reales
# predicción 0 1
# 0 305 201
# 1 301 443
#Eficiencia de la clasificación 59.8%

La gráfica que se obtiene del SVM con kernel lineal es:

SVM_Mezcla3_ker-lineal

Haciendo la comparación entre datos originales y las predicciones se tiene lo siguiente:

Orig3_vs_Ker-lineal

Se aprecia que no es muy favorable usar el kernel lineal con estos datos, la eficiencia de la clasificación es del 59.4%

Para comparar con otro kernel uso el radial para hacer la clasificación.

#Kernel Radial

set.seed(1)
tune.out=tune(svm,factor(y)~.,data=X,kernel="radial",sacale=FALSE,type='C-classification',ranges=list(cost=c(0.001,0.01, 0.1, 1,5,10,100)))
tune.out
summary(tune.out)
bestmod=tune.out$best.model
summary(bestmod)
#Informa del tipo de modelo, del que puede considerarse como el mejor modelo
Pred=predict(bestmod,X)
table(predicción=Pred,Valores_reales=X$y)
# Valores_reales
# predicción 0 1
# 0 603 7
# 1 3 637
#Eficiencia de la clasificación 99.2%, con kernel radial

La gráfica que se obtiene con este kernel es la siguiente:

SVM_Mezcla3_ker-radial

Esta imagen muestra que las curvas que se obtienen con el kernel radial son parecidas a las de los datos originales. Haciendo la comparación entre los datos y las estimaciones se obtiene lo siguiente:

Orig3_vs_Ker-Radial

Se aprecia que es muy similar y más aún la eficiencia es del 99.2%, lo cual comparado con el lineal es sumamente mejor usar el kernel radial.

Ejemplo con datos de correos clasificados

Los datos se encuentran cargados en la librería kernelabs, con la cual también se puede estimar la máquina de soporte vectorial. Las comparaciones las hago con varias técnicas:

  • Naive Bayes
  • LDA y Regresión Logística
  • Redes Neuronales
  • Árboles de decisión y random forest

Antes de eso aplico con los cuatro kernel disponibles en svm de e1071 para elegir el que mejor clasificación realiza en la muestra de prueba (test set).

Primero reviso los datos.

#Exploración breve de los datos
#Se cargan primero las librerías

library(kernlab)
library(e1071)
library(MASS)
library(nnet)
librery(randomForest)
data(spam)

#Datos
head(spam)
dim(spam)
#4601 58
str(spam)

#############################################
#Preparación de datos

index=sort(sample.int(4601,1150))
spam.train=spam[-index,]
spam.test=spam[index,]

Se puede hacer algo más al explorar la información, diseñar algunas gráficas o usar alguna técnica de estadística de exploración.

Lo que se hace en el código anterior es dividir los datos entre una conjunto de entrenamiento y un conjunto de prueba.

Aplico los cuatro kernel para elegir el que mejor clasificación realiza.

#SVM con diferentes kernel
########################################
#Kernel lineal

svmfit=svm(type~.,data=spam.train,kernel="linear",cost=0.1,type='C-classification' )
summary(svmfit)
Pred=predict(svmfit,spam.test)
table(predicción=Pred,Valores_reales=spam.test$type)
#         Valores_reales
#predicción nonspam spam
#   nonspam 646 40
#   spam     40 424
#Eficiencia del modelo 93.04%
#######################################
#Kernel radial
svmfit=svm(type~.,data=spam.train,kernel="radial",cost=0.1,type='C-classification' )
summary(svmfit)

set.seed(1)
tune.out=tune(svm,type~.,data=spam.train,kernel="radial",type='C-classification',ranges=list(cost=c(0.001,0.01, 0.1, 1,10,100)))
tune.out
summary(tune.out)
bestmod_r=tune.out$best.model
summary(bestmod_r)
Pred=predict(bestmod_r,spam.test)
table(predicción=Pred,Valores_reales=spam.test$type)
#           Valores_reales
#predicción nonspam spam
#   nonspam     652 36
#   spam         34 428
#Eficiencia del modelo 93.91%
#################################################
#Kernel sigmoid
svmfit=svm(type~.,data=spam.train,kernel="sigmoid",cost=0.1,type='C-classification' )
summary(svmfit)
set.seed(1)
tune.out=tune(svm,type~.,data=spam.train,kernel="sigmoid",type='C-classification',ranges=list(cost=c(0.001,0.01, 0.1, 1,10,100)))
tune.out
summary(tune.out)
bestmod_s=tune.out$best.model
summary(bestmod_s)
Pred=predict(bestmod_s,spam.test)
table(predicción=Pred,Valores_reales=spam.test$type)
#           Valores_reales
#predicción nonspam spam
#    nonspam    649 78
#    spam        37 386
#Eficiencia del modelo 90%
################################################
#Kernel polynomial
svmfit=svm(type~.,data=spam.train,kernel="polynomial",cost=0.1,type='C-classification' )
summary(svmfit)
set.seed(1)
tune.out=tune(svm,type~.,data=spam.train,kernel="polynomial",type='C-classification',ranges=list(cost=c(0.001,0.01, 0.1, 1,10,100)))
tune.out
summary(tune.out)
bestmod_p=tune.out$best.model
summary(bestmod_p)
Pred=predict(bestmod_p,spam.test)
table(predicción=Pred,Valores_reales=spam.test$type)
#           Valores_reales
#predicción nonspam spam
#   nonspam     661 53
#   spam         25 411
#Eficiencia del modelo 93.2%

Se observa en las tasas de error que el mejor kernel de svm es el radial, así que considero este contra el cual comparo las otras técnicas.

Si bien no explico a detalle lo que hace el código en breve es estimar vía tune() la elección del mejor valor de cross-validation de varios valores.

Una técnica clásica para clasificar correos es Naive Bayes, la cual considera ciertas propiedades teóricas que permiten hacer un calculo rápido. Uso NaiveBayes de la librería e1071.

#Naive Bayes vs SVM elegido
NB=naiveBayes(type~.,data=spam.train)
Pred=predict(NB,spam.test)

table(predicción=Pred,Valores_reales=spam.test$type)
# Valores_reales
#predicción nonspam spam
# nonspam 362 18
# spam 324 446
#Se tiene una eficiencia del 70.2% de bien clasificados

Resulta que es mejor el modelo SVM ante el Naive Bayes por un porcentaje altamente mejor.

#Regresión Logística-LDA vs SVM elegido

Mod.lm=glm(type~.,data=spam.train,family='binomial')
summary(Mod.lm)
Mod.lmprob=predict(Mod.lm,newdata=spam.test,type="response")
Mod.pred=rep('nonspam',1150)
Mod.pred[Mod.lmprob>0.5]='spam'
table(Predicción=Mod.pred,Valores_Reales=spam.test$type)
#           Valores_Reales
#Predicción nonspam spam
#   nonspam 643 39
#   spam     43 425
#Predicción 92.86%
############################################
#Método LDA
lda.fit=lda(type~.,data=spam.train)
lda.pred=predict(lda.fit,spam.test)
lda.class=lda.pred$class
table(Predicción=lda.class,Valores_Reales=spam.test$type)
#lda.class nonspam spam
#   nonspam    650 91
#   spam        36 373
#Tasa de eficiencia 88.9%

Se tienen que el porcentaje de eficiencia en la clasificación de SVM 93.91% es mayor de casi 1% mejor que los métodos lineales, en este caso la regresión logística resulta ser la mejor entre las técnicas lineales.

Sin explicar detalles elijo un modelo de redes neuronales para clasificar, la librería es nnet y la función tienen el mismo nombre.

#Redes Neuronales
red=nnet(type~.,data=spam.train,size=2)
Pred=predict(red,spam.test,type="class")

table(Predicción=Pred,Valores_Reales=spam.test$type)
# Valores_Reales
#Predicción nonspam spam
# nonspam 650 29
# spam 36 435
#Tasa de eficiencia 94.34%

Resulta que el modelo de redes neuronales resulta ser mejor casi 0.5% mejor que el modelo de SVM. Lo cual faltaría hacer una revisión para definir el mejor modelo de redes neuronales y quizás la mejora es mejor o nula pero tendrías los valores óptimos de los parámetros.

#Dos técnicas de arboles
# La primera es por árboles y se hace cross-validation para elegir el mejor modelo
 
tree.spam=tree(type~.,spam.train)
set.seed (3)
cv =cv.tree(tree.spam,FUN=prune.misclass )

#Gráfica de Cross-validation para elegir el mejor valor de del parámetro best
par(mfrow =c(1,2))
plot(cv$size,cv$dev,type="b",col="2")
plot(cv$k,cv$dev,type="b",col="3")

prune.tree=prune.misclass(tree.spam,best=8)
tree.pred=predict(prune.tree,spam.test,type="class")
table(tree.pred,Valores_Reales=spam.test$type)

#         Valores_Reales
#tree.pred nonspam spam
#  nonspam 644 58
#  spam     42 406
#Tasa de clasificación 91.3%
#################################################
#Modelo de Random Forest

set.seed (1)
rf.spam =randomForest(type~.,data=spam.train,mtry=7, importance =TRUE)
rf.pred_spam=predict(rf.spam,newdata=spam.test,type="class")

table(rf.pred_spam,Valores_Reales=spam.test$type)
#            Valores_Reales
#rf.pred_spam nonspam spam
#     nonspam 663 29
#     spam     23 435
#Tasa de clasificación 95.4%

Revisando la tasa de clasificaciones correcta por medio de las técnicas de árboles de decisión , resulta ser la mejor técnica para clasificar los datos de los correos. Se obtienen que la técnica de Random Forest resulta clasificar correctamente el 95.4%, comparado con la red neuronal y SVM es mejor casi en un 2%.

La única intensión de este ejemplo es hacer una comparación entre diversas técnicas de clasificación, de cierto modo hacer un ejemplo sencillo y que se puede replicar. Ha esto le faltaría agregar algunos gráficos que muestran los resultados o hacer algunas variaciones de los parámetros requeridos en cada algoritmos para determinar con mejor eficiencia el mejor modelo.

Espero ilustre de manera breve la variedad de técnicas y permita visualizar que al final se requiere probar con muchos modelos para elegir uno entre todos los realizados.