por ·Sin comentarios

El 9 de octubre nos publicaron un paper en Intelligent Data Analysis, colaboración con Mimi Chong y Matt Davison de la Universidad de Western Ontario, en Canadá. El paper se titula «How much effort should be spent to detect fraudulent applications when engaged in classifier-based lending?» y analizamos teóricamente cuál es el costo real de las mentiras en las aplicaciones crediticias. El abstract sale más abajo.

El link a la revista es este.

Abstract: Credit scoring is an automated, objective and consistent tool which helps lenders to provide quick loan decisions. It can replace some of the more mechanical work done by experienced loan officers whose decisions are intuitive but potentially subject to bias. Prospective borrowers may have a strong motivation to fraudulently falsify one or more of the attributes they report on their application form. Applicants learn about the characteristics that are used to build credit scoring models, and may alter the answers on their application form to improve their chance of loan approval. Few automated credit scoring models have considered falsified information from borrowers. We will show that sometimes it is profitable for financial institutions to spend money and effort to identify dishonest customers. We will also find the optimal effort that banks should spend on identifying these liars. Furthermore, we will show that it is possible for liars to eventually adjust their lies to escape from credit checks. The proposed issue will be studied using simulated data and discriminant analysis. This research can help lending financial institutions to reduce risk and maximize profit, and it also shows that it is feasible for customers to lie intelligently so as to evade credit checks and get loans.

por ·Sin comentarios

En este post mostraré algunos trucos que he ido aprendiendo al utilizar Random Forests (RF), una de mis técnicas de clasificación favoritas, para mejorar el poder de clasificación que poseen. Este método, desarrollado por Leo Breiman en el año 2001 tiene sus orígenes en esta publicación en Machine Learning y se ha convertido, especialmente en los últimos cinco o seis años, en una de las preferidas por científicos y practitioners. La técnica es tremendamente poderosa, es muy fácil de usar (en comparación a las redes neuronales por ejemplo), y entrega bastante información sobre los resultados al usuario.

La lógica de los RF es generar muchos árboles de decisión sobre los datos de entrada, utilizando una cantidad reducida de las variables que se encuentran disponibles. La técnica explora así diferentes subsegmentos del espacio de entrada, y aprende patrones muy complejos aprovechando la diversidad de los clasificadores que entrena. Una prueba de su éxito se ha visto en las competencias de Kaggle: es una de las dos técnicas más importantes que identifican los administradores del sitio.

En este ejemplo entrenaré un RF para un dataset del repositorio de la UCI, y luego mostraré algunos trucos: cómo acelerar el entrenamiento utilizando programación paralela para aprovechar mejor los recursos disponibles y cómo graficar una curva ROC utilizando

ggplot2

. En un post futuro mostraré cómo balancear las muestras para solucionar los problemas de desbalance y finalmente cómo graficar la importancia de las variables y analizar estos resultados.

Entrenamiento Secuencial

Primero importemos los datos. Para este ejemplo usaré la base de datos de Bank Marketing, de Moro et. al. Esta base de datos posee 41.188 instancias, 20 variables predictoras, y tiene como objetivo predecir si el cliente toma o no un depósito a plazo. Usaré R para descargar directamente el archivo, descomprimirlo, y cargar los datos. Al final agrego una variable adicional para los test (muestra holdout).

# Descargar el archivo. Se debe usar "method=curl" si la página usa https.
url <- "https://archive.ics.uci.edu/ml/machine-learning-databases/00222/bank-additional.zip";
download.file(url,"bank-additional.zip", mode="wb",method="curl");
# Descomprimir el archivo y cargar la tabla.
unzip("bank-additional.zip");
bank.data <- read.table("bank-additional-full.csv", sep=";", header = TRUE);
bank.data$if.test <- rbinom(nrow(bank.data), 1, 0.7) # Agrega variable para test

Con esto el archivo está cargado en la variable

bank.data

y ya se puede utilizar. Entrenemos un primer Random Forest y veamos sus resultados. Para entrenar los modelos está el paquete

randomForest

de Breiman mismo, y además utilizaré los paquetes

verification

,

rbenchmark

,

e1071

y

caret

para analizar los resultados. Carguemos los paquetes o instalémoslos si es necesario.

# library(..., logical.return = TRUE) devuelve falso si el paquete no existe.
if(!library(randomForest, logical.return = TRUE)) {
  install.packages('randomForest');
}
if(!library(verification, logical.return = TRUE)) {
  install.packages('verification');
}
if(!library(caret, logical.return = TRUE)) {
  install.packages('caret');
}
if(!library(e1071, logical.return = TRUE)) {
  install.packages('e1071');
  library(e1071)
}
if(!library(rbenchmark, logical.return = TRUE)) {
  install.packages('rbenchmark');
  library(rbenchmark)
}

Ahora podemos correr un primer Random Forest.

# Crea fórmula con todas las variables. Con los paste evito escribir de más.
vars.bank <- colnames(bank.data)
formula.rf <- as.formula(paste0("y ~", paste(vars.bank[1:20], collapse = "+")))

# Función para facilitar futuras ejecuciones.
rfNaive <- function(x){randomForest(formula.rf, data = x, 
                                     subset = (if.test == 0), ntree = 500,
                                     importance = TRUE)
}
# Entrena sin ninguna precaución. Corre usando 650MB en RAM y en
# 36 segundos en un procesador core i7-3770k.
rf.naive <- rfNaive(bank.data)

# Calcula AUC y Crea matriz  de confusión con corte en 0.5.
prob.rfNaive <- predict(rf.naive, newdata = bank.data[bank.data$if.test == 1 , ],
                        type = "prob")[, 2]
roc.area(as.numeric(bank.data$y[bank.data$if.test == 1] == "yes"), prob.rfNaive)
confusionMatrix(round(prob.rfNaive),  as.numeric(bank.data$y[bank.data$if.test == 1] == "yes"))
# AUC = 0.9426409
# Matriz de confusión:
#               Reference
# Prediction     0     1
#         0   24711  1546
#         1   939  1703

Entrenando en Paralelo

Los Random Forest son altamente paralelizables pues cada árbol puede entrenarse por si solo, y al final del entrenamiento basta con juntar cada uno de los árboles que se entrenaron. El paquete

randomForest

toma ventaja de eso y entrega directamente una manera de combinar resultados cuando se entrenan bosques en paralelo.

Para paralelizar en R es necesario registrar un backend que es el que se encarga de comunicarse con los procesadores, y además usar un paquete que permita paralelizar el código. Hay distintos paquetes con distintos niveles de profundidad en cuánto a control (y cuanta dificultad) entregan al usuario. Lo que me ha resultado más fácil es:

  • Usar como backend ya sea el paquete doMC (Linux) o el paquete doParallel(Windows). Lamentablemente no existe (que yo sepa) un paquete que funcione en ambos sistemas.
  • Usar el paquete foreach para realizar la paralelización.

Las opciones de estos paquetes son amplias, y escapan el alcance de este post, así que me enfocaré en lo necesario para correr randomForests. Para todos los detalles de paralelizar R se puede comenzar por la viñeta de foreach.

¿Qué es paralelizar? La gran mayoría de los computadores modernos traen más de un procesador (núcleo) incorporado para correr procesos. La ventaja de tener más de un procesador radica en que el usuario usa múltiples programas y servicios en su equipo en cualquier momento, y el tener más de un núcleo permite poder realizar varias actividades de forma simultánea. Nosotros podemos tomar ventaja de esto para generar códigos más eficientes, que utilicen todo el poder computacional que tenemos. Para paralelizar un Random Forest la estrategia a seguir será simplemente decirle al programa que entrene una cierta cantidad de árboles en cada núcleo, y después combine los resultados en un único bosque aleatorio.

La computación paralela tiene una propiedad que es importante tener presente: dependiendo de cuán eficiente es la implementación del código, puede multiplicar la cantidad de memoria utilizada por la cantidad de procesadores que se utilizan. Así, si una aplicación utiliza 500MB en RAM en un núcleo, perfectamente puede utilizar 2GB si se utilizan cuatro, por ejemplo. Por lo mismo, mucho cuidado al correr este ejemplo, ajusten los valores de los parámetros acorde a la memoria RAM que tengan disponibles.

Para paralelizar nuestro Random Forest, lo primero es determinar cuántos núcleos poseemos y cuántos árboles queremos entrenar. En este ejemplo estamos entrenando 500 árboles, y en mi PC tengo ocho núcleos y suficiente RAM, así que los usaré todos. El número de procesadores se puede obtener googleando el tipo de procesador que tienen, o en Windows con el backend parallel pueden correr

detectCores()

en R, o en Linux el backend doMC utiliza la mitad de los núcleos disponibles de forma predeterminada.

Con eso listo ahora podemos paralelizar. Le diremos a R que entrene sub-bosques, y después los combine. El código es el siguiente:

cores <- 4 # Depende de cada computador.
ntrees <- 500 # Decisión del usuario.
library(doMC)
registerDoMC(cores) # Registra el backend. Solo en Linux.
# En Windows:
# library(doParallel)
# registerDoParallel(cores = cores)

rfParallel <- function(x){
  foreach(ntree.iter = rep(ceiling(ntrees / cores), cores), 
          .combine = combine, .packages = "randomForest") %dopar% {
            randomForest(formula.rf, data = x, ntree = ntree.iter,
                         do.trace = FALSE, subset = (if.test == 0),
                         keep.forest = TRUE, importance = TRUE)
          }
}

rf.parallel <- rfParallel(bank.data)

# Calcular medidas efectividad
prob.rfParallel <- predict(rf.parallel, newdata = bank.data[bank.data$if.test == 1 , ],
                        type = "prob")[, 2]
roc.area(as.numeric(bank.data$y[bank.data$if.test == 1] == "yes"), prob.rfParallel)
confusionMatrix(round(prob.rfParallel),  
                as.numeric(bank.data$y[bank.data$if.test == 1] == "yes"))
# AUC = 0.9439
# Matriz de confusión:
#               Reference
# Prediction     0     1
#       0      24665  1584
#       1        905  1673
# Accuracy: 91.37%
# Balanced Accuracy: 73.91%

Como podemos ver, los resultados son equivalentes entre los dos modelos. Existe una pequeña caída en la capacidad discriminante, pero esto es sólo por razones aleatorias.

Los parámetros de la función

foreach

son los importantes:

  •  ntree.iter = rep(ceiling(ntrees / cores), cores)

    : Indica que cada iteración de foreach entrene un total de 500 divido por la cantidad de núcleos que disponen.

  • .combine = combine

    : Indica a 

    foreach

    que utilice la función «combine» del paquete 

    randomForest

    para combinar los resultados. Esta función permite unir bosques aleatorios en bosques más grandes.

  • .packages='randomForest'

    : Indica a 

    foreach

    que, en cada núcleo donde correrá la función, cargue el paquete 

    randomForest

    .

  • %dopar%

    : Indica a 

    foreach

    que corra en paralelo utilizando el backend que se encuentre registrado (

    doMC

    en este caso).

Ahora comparemos el tiempo de corrida. Repetiré tres veces el entrenamiento utilizando el excelente paquete

rBenchmark

.

benchmark(rfNaive(bank.data), rfParallel(bank.data), 
          order = "elapsed", replications = 3,
          columns = c("test", "replications", "elapsed", "relative",
                      "user.self"))

#         test           replications elapsed relative user.self
#2 rfParallel(bank.data)       3       77.760   1.00     2.204
#1    rfNaive(bank.data)       3      148.557   1.91   138.380

Este ejemplo lo corrí en mi laptop con dos núcleos reales y cuatro virtuales. La programación en paralelo corre casi dos veces más rápido que el código secuencial. Esto es en general lo que ocurre: existe una ganacia que es casi lineal con el número de núcleos físicos que se disponen.

 

Consejo 2: Graficar las curvas ROC con ggplot2

Para finalizar este post mostraré algo simple para quienes reportan los resultados de modelos. Las curvas ROC de R en general se ven bastante mal, incluso en el mismo paquete

verification

indican que existen mejores alternativas para graficar.

El paquete más potente para gráficos en R corresponde a ggplot2, pero utilizarlos es bastante complejo. El paquete utiliza el concepto de estéticas, o variables que se asignan a distintas partes del gráfico. Toma un poco de tiempo acostumbrarse a esto, pero una vez que se logra los resultados son muy buenos.

Para graficar una curva ROC, es necesario primero que nada de disponer un data.frame tal que muestre los puntos de la curva ROC y a quién se refiere. Luego creamos el gráfico, asignamos la variable que identifica al modelo tanto a la estética de grupos como a la de colores, y luego le damos los valores correspondientes. El código es el siguiente:

# Crea estructura con curva ROC
roc.rfParallel <- roc.plot(as.numeric(bank.data$y[bank.data$if.test == 1] == "yes"), 
                           prob.rfParallel)
roc.rfParallel <- as.data.frame(roc.rfParallel$plot.data)
roc.rfNaive <- roc.plot(as.numeric(bank.data$y[bank.data$if.test == 1] == "yes"), 
                           prob.rfNaive)
roc.rfNaive <- as.data.frame(roc.rfNaive$plot.data)

temp <- data.frame(fpr = roc.rfParallel$V3, var = 'En Paralelo', value = roc.rfParallel$V2)
temp <- rbind(temp, data.frame(fpr = roc.rfNaive$V3, var = 'Secuencial', value = roc.rfNaive$V2))

# Leyenda del gráfico
legend.roc <- rep("",2)
legend.roc[1] <- paste0("En Paralelo: AUC = ", round(with(bank.data, roc.area(as.numeric(y[if.test == 1] == "yes"),
                                                                              prob.rfParallel))$A,3))
legend.roc[2] <- paste0("Secuencial: AUC = ", round(with(bank.data, roc.area(as.numeric(y[if.test == 1] == "yes"),
                                                                                 prob.rfNaive))$A,3))

breaks <- c("En Paralelo", "Secuencial")

# Gráfico
library(ggplot2)
library(grid)

p <- (ggplot(temp, aes(fpr, value, group = var, color = var)) 
      + geom_line(size = 0.5, aes(linetype = var))+ theme_bw()
      + scale_linetype_manual(name = 'leyenda', values=c("dotted", "solid"), 
                              labels = legend.roc, breaks = breaks)
      + scale_color_discrete(name = 'leyenda', breaks = breaks, labels = legend.roc)
      + scale_x_continuous("False Positive Rate (1-Specificity)")
      + scale_y_continuous("True Positive Rate (Sensitivity)")
      + coord_fixed()
      + theme(legend.justification=c(1,0), legend.position=c(1,0),
              legend.title = element_blank(), 
              legend.text = element_text(size = 14),
              legend.key.width = unit(3, "cm"),
              legend.key = element_blank()
      )
)
plot(p)

El resultado es:

 

Curva ROC Random Forests

Curva ROC Random Forests

Como se observa, ambos modelos son idénticos. Esto puede cambiar en una comparación entre distintos tipos de modelos.

Este post mostró tanto como entrenar los bosques aleatorios en paralelo y como generar curvas ROC con

ggplot2

. Cualquier comentario, por twitter o directamente acá en el blog. ¡Gracias por leer!

Actualización: Una versión anterior de este post tenía algunos paquetes faltantes en el código. ¡Gracias Pablo!

por ·Sin comentarios

En los últimos meses he estado tratando de aprender a optimizar mi código en R. En general, R funciona bastante bien, pero hay varias aplicaciones donde se puede poner muy, pero muy lento, en particular los loops. Un búsqueda rápida muestra que existen muchos artículos que discuten el tema con bastante profundidad.

Uno de los temas que encontré y llamaron mi atención fue la capacidad de integrar R con otros lenguajes de programación, en particular con algunos de más bajo nivel, como C, Fortran o C++. La integración con C y Fortran está ya solucionada directamente en R, a través de la función .Call, pero esta forma de llamar a código externo tiene la gran desventaja que todas las funciones deben ser declaradas como void, es decir, debemos modificar los valores entregados usando directamente punteros a los objetos en memoria. Esto es… engorroso, al menos lo fue en mis primeras experiencias programando.

La solución: usar C++. Es un lenguaje moderno, que corre eficientemente los loops, y que tiene muchas herramientas que facilitan la vida al programador. Para integrar R y C++ existe el excelente paquete Rcpp, el cual permite crear funciones de forma dinámica y sencilla en R, y que además permite usar objetos inherentes a R (vectores numéricos, vectores de caractéres, dataframes, etc.) en C++, pero aún este paquete tiene una falencia: no soporta operaciones matriciales eficientes. Para ello, un plugin del mismo paquete llamado RcppArmadillo ofrece una interfaz entre Rcpp y la librería Armadillo que entrega herramientas para operaciones matriciales en C++. Con esto tenemos todo lo que necesitamos.

En este artículo programaré un código sencillo para obtener un clúster de K-Medias y compararé las implementaciones en código en R, la implementación oficial (en C) y una implementación propia en C++. El algoritmo es bastante sencillo: inicializa vectores de forma aleatoria, para luego ir ajustando los clusters en cada iteración hasta alcanzar convergencia.

Implementación (ineficiente) en R.

kmeans.R <- function(x,k,eps = 10e-6){
  # k es la cantidad de clusters, x es la matriz de datos.
  # eps es la tolerancia máxima entre una iteración y otra.
  require(pdist)
  cluster <- rep(0,nrow(x)) # Cluster actual del caso.
  # Inicializa centroides.
  centroids <- matrix(runif(nrow(x) * k), nrow = k, ncol = ncol(x))
  # Guarda última iteración
  centroids.old <- matrix(0, nrow = k, ncol = ncol(x))
  eps.vec <- rep(eps,k)
  dists.centroids <- 1000
  while(any(dists.centroids > eps.vec)){
    centroids.old <- centroids
    # Matriz con las distancias entre casos y vectores.
    dists.x <- as.matrix(pdist(x, centroids)) 
    min.x <- apply(dists.x, 1, min)
    # Obtiene indice del menor valor, i.e., el cluster.
    cluster <- apply(cbind(dists.x,min.x), 1, 
                     function(x) which(x[1:k] == x[k+1])[1])     
    for(i in 1:k){
      # Calcula el centroide.
      centroids[i,] <- apply(x[cluster == i, ], 2, mean) 
    }
    dists.centroids <- diag(as.matrix(pdist(centroids, centroids.old)))
  }
  # Devuelve cluster y centroides.
  return(list(cluster = cluster, centroids = centroids)) 
}

Para probar el código, en R simularemos una matriz sencilla con dos distribuciones. La primera una normal de estándar y la segunda una normal de media 5 y desviación estándar 1. Posteriormente aplicamos la función y graficamos. Noten que estoy usando el paquete ‘pdist’, solo para facilitarme la vida cuando calculo las distancias entre matrices. Este paquete tira algunas warnings en la última iteración, simplemente para avisar que convergió, así que ignórenlas.

library(pdist)
x <- rbind(matrix(rnorm(1000 * 2, mean = 0), nrow = 1000, ncol = 2), 
           matrix(rnorm(1000 * 2, mean = 5), nrow = 1000, ncol = 2))
x.kmeans.R <- kmeans.R(x, 2)
x.kmeans.R$centroids
# Resultado:
#            [,1]       [,2]
# [1,] 0.04097412 0.04618219
# [2,] 5.02960813 4.97157822
plot(x, col=ifelse(x.kmeans.R$cluster == 1, "red", "green"))

El resultado es el siguiente:

Datos en dos segmentos.

Implementación en C++.
Implementemos ahora la misma función en C++. Para ello usaremos Rcpp y RcppArmadillo. Primero, es necesario instalar los paquetes:

install.packages(c('Rcpp', 'RcppArmadillo'))

.

Ahora, existen varias rutas para implementar la función. Es posible escribir directamente la función en un gran string de texto y compilarla en línea, o se puede guardar un archivo anexo que tenga la función programada. Seguiremos este camino. En un nuevo archivo (kmeansCpp.cpp) programamos la función de kmedias.

#include <RcppArmadillo.h>
using namespace Rcpp;

// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
List kmeansCpp(NumericMatrix xa, int k, double eps = 10e-3) {
  arma::mat x(xa.begin(), xa.nrow(), xa.ncol(), false);
  // Inicializa centroides de forma uniforme en [0-1]
  arma::mat centroids(k, xa.ncol(), arma::fill::randu);
  // Incializa con puros 0's.
  arma::mat centroids_old(k, xa.ncol(), arma::fill::zeros);
  // Crea vector fila tamaño 1xk.
  arma::rowvec eps_vec = eps * arma::ones<arma::rowvec>(k);
  arma::rowvec dists_centroids = 100 * arma::ones<arma::rowvec>(k);
  // Crea matriz tamaño casos x k.
  arma::mat dists_x(xa.nrow(), k);
  LogicalMatrix temp(xa.nrow(), k);
  arma::mat cluster(xa.nrow(), k);
  while(arma::any(dists_centroids > eps_vec)){
    centroids_old = centroids;
    // Calcula matriz de distancias
    for(int i = 0; i<k; i++){
      dists_x.col(i) = sqrt(sum((x - (arma::ones(xa.nrow(), 1) 
                                  * centroids.row(i))) % 
                      (x - (arma::ones(xa.nrow(), 1) * centroids.row(i))), 1));
    }
    // Obtiene indice del menor valor, i.e., el cluster.
    arma::colvec min_x = min(dists_x, 1);
    temp = wrap((min_x * arma::ones(1, k) == dists_x));
    cluster = as<arma::mat>(temp);
    // Calcula el centroide.
    for(int ia = 0; ia < k; ia++){
      centroids.row(ia) = mean(x.rows(find(cluster.col(ia) == 1)), 0);
    }
    dists_centroids = (sqrt(sum((centroids - centroids_old) % 
                      (centroids - centroids_old), 1))).t();
  }
  // Devuelve cluster y centroides.
  return List::create(Named("cluster") = wrap(cluster), 
                      Named("centroids") = wrap(centroids));
}

Analicemos un poco el programa:

  • La línea #include <RcppArmadillo.h> le indica al compilador de C++ que importe la cabecera que nos permitirá usar todas las herramientas para la interfaz con R.
  • using namespaces Rcpp nos permite utilizar todas las funciones de Rcpp sin tener que colocar Rcpp::NombreFuncion cada vez que la llamamos. También puede ser útil colocar using namespaces arma para ahorrarse los arma::.
  • // [[Rcpp::depends(RcppArmadillo)]] y // [[Rcpp::export]] le dicen a Rcpp que tiene que ocupar Armadillo y que debe exportar la función a R.

Luego viene la función en si. Tanto Rcpp como Armadillo poseen clases que permiten manejar vectores, matrices, listas, etc. En este caso, es importante notar que creamos una función que recibe y entrega objetos de Rcpp. La función recibe una matriz numérica (NumericMatrix, objeto de Rcpp), un entero y un double y entrega una lista de R. Los objetos de Rcpp son, a grandes rasgos, matrices (Matrix), vectores (Vector), listas (List) y Dataframes (Dataframe); mientras que los dos primeros pueden ser numéricos (Numeric), enteros (Integer), literales (String) o booleanos (Logical). Así, una matriz numérica corresponde a un objeto NumericMatrix, mientras que un vector booleano será un LogicalVector. Cada objeto tiene diversos métodos que se le pueden aplicar, pero en general las operaciones matemáticas funcionan correctamente. Para una descripción detallada, este tutorial es un buen punto de partida.

Dentro de la función, las primeras líneas inicializan los objetos que usaremos en la iteración. Por ejemplo, arma::mat x(xa.begin(), xa.nrow(), xa.ncol(), false); indica que crearemos una matriz x, objeto de Armadillo (los objetos más usados son vec, para vector numérico y mat para matrices numéricas) reutilizando la memoria que apunta a la matriz numérica de Rcpp xa, nuestro input. Los demás objetos se inicializan a partir de los tipos que usaremos.

El siguiente paso es el loop principal, que sigue la misma estructura que el loop de R, pero ahora utilizando las funciones propias de Armadillo o Rcpp, dependiendo de cuáles sean los objetos sobre los que se las estamos aplicando. Por ejemplo, el operador % corresponde a la multiplicación término a término entre matrices de Armadillo. Finalmente, devolvemos una lista de R con el cluster (ahora como una matriz con índices) y los centroides. La función wrap() permite transformar objetos de Armadillo en objetos de Rcpp fácilmente.

Para poder utilizar la función en R en nuestra sesión corremos el comando Rcpp::sourceCpp("kmeansCpp.cpp"). Con ello la función kmeansCpp() está disponible para su uso.

Comparando resultados

Comparemos los códigos, para ver cuánto ganamos por traspasar la función a C++. Utilizaré el paquete rbenchmark que permite rápidamente comparar distintas funciones en cuanto a su tiempo de ejecución. El código de comparación completo queda entonces:

library(RcppArmadillo)
library(rbenchmark)
Rcpp::sourceCpp("kmeansCpp.cpp")
x <- rbind(matrix(rnorm(1000 * 2, mean = 0), nrow = 1000, ncol = 2), 
           matrix(rnorm(1000 * 2, mean = 5), nrow = 1000, ncol = 2))
benchmark(kmeans(x,2), kmeansCpp(x,2), kmeans.R(x,2), order = "elapsed",
          columns = c("test", "replications", "elapsed", "relative", 
                      "user.self"))
# Resultados:
#              test replications elapsed relative user.self
# 2 kmeansCpp(x, 2)          100   0.040    1.000     0.036
# 1    kmeans(x, 2)          100   0.059    1.475     0.056
# 3  kmeans.R(x, 2)          100   7.775  194.375     7.752

Nuestra implementación en C++ corre en 0.036 segundos, 194 veces más rápido (!!!) que la implementación en R puro y corre además casi un 50% más rápido que la implementación en C que trae incorporado R. El resultado habla por si solo, vale la pena traspasar códigos con muchas iteraciones a C++, pues es más eficiente tanto en uso de memoria como en velocidad. La diferencia con kmeans de R es engañosa, pues no realizamos ningún control de errores ni depuraciones a las entradas que probablemente sí realiza esta función, por lo que las ganancias se deben a eso probablemente.

Así que ya saben, si quieren que su código vuele, bajen un poco de nivel y utilicen código en C++. R ofrece todas las herramientas para lograrlo de forma sencilla.

por ·2 comentarios

Primero que nada, ¡bienvenidos a mi primer post en el blog de mi página personal! Instalé WordPress en un subdirectorio de mi servidor principal, por lo mismo no se ve tan bonito como la página original. Con el tiempo iré puliendo cómo se va viendo 🙂

Ahora, a lo que vinimos. Hace un par de fin de semanas, e inspirando en esta noticia, me nació el interés de crear mapas dinámicos que pudieran, usando código, visualizar Chile. Más de alguna vez he tenido ganas de mostrar algo que pase en el país a través de mapas, pero la sola idea de lanzarme a utilizar Illustrator u otra herramienta para diseñadores me daba náuseas. Me gusta programar, no dibujar, así que ¡¿cómo uso código dammit?!

Partamos. Primero que nada, los créditos. Este post está basado casi enteramente en el tutorial escrito, para Alemania, por Mark Heckman. Pueden ir a leer sus instrucciones allí si lo desean, yo adapté todo para el caso chileno.

¿Qué es un mapa? Para nosotros un mapa es simplemente una colección de coordenadas que, graficadas en un plano, representan una región geográfica real. Nosotros buscamos un mapa político, pues queremos la división regional.

Para el caso chileno tenemos dos grandes recursos: los Mapas Vectoriales de la Biblioteca del Congreso Nacional, en formato shapefile y utilizados en la parte dos de este tutorial, y el proyecto Global Administrative Areas. Este último tiene como objetivo obtener todos los mapas políticos del mundo al mayor de nivel de detalle posible en diversos formatos y están disponible los chilenos inmediatamente en R. En el caso que quieran mapas geográficos, las Generic Mapping Tools entregan recursos excelentes, pero escapa a los alcances de este post.

Ahora a generar el mapa. Para este ejemplo generaremos un mapa sencillo con la tasa de default en créditos por región, y asumiré que se manejan con el software R, que a mi gusto entrega la mayor cantidad de recursos estadísticos de forma gratuita.

Paso 1: Instalar los paquetes necesarios y cargarlos.

install.packages('sp') #Permite leer los mapas.
install.packages('RColorBrewer') #Paletas de colores para mapas.
load(sp)
load(RColorBrewer)

Paso 2: Descargar el mapa de Chile: 

con <- url('http://biogeo.ucdavis.edu/data/gadm2.8/rds/CHL_adm1.rds')
print(load(con))
close(con)

Este archivo contiene los datos de Chile a nivel de provincia actualizados al año 2007, por lo tanto figuran solamente 13 regiones.

Paso 3: Corregir las regiones faltantes: El objeto «gadm» posee una serie de variables que describen al país. NAME_1 incluye el nombre de las regiones, y NAME_2 el de las provincias. Para obtener las 15 regiones simplemente debemos asignar las provincias de Arica y Parinacota a la región XV y la provincia de Valdivia a la región de Los Ríos. Por simpleza esto lo realicé a mano. Una vez arreglado, debemos transformar a factores las dos variables.

gadm$NAME_1[gadm$NAME_2 == "Valdivia"] <- "Los Ríos"
gadm$NAME_1[gadm$NAME_2 == "Arica"] <- "Arica y Parinacota"
gadm$NAME_1[gadm$NAME_2 == "Parinacota"] <- "Arica y Parinacota"
gadm$NAME_1 <- as.factor(gadm$NAME_1)
gadm$NAME_2 <- as.factor(gadm$NAME_2)

Paso 4: Cargar los datos a corregir

Como estamos usando los datos por provincia, necesitamos una variable adicional (que puede estar dentro o no de gadm) tal que a cada elemento de NAME_1 (las regiones) le asocie un valor que será el que graficaremos.  Este archivo posee los datos por default por región, extraído a partir de las colocaciones por región e institución a Marzo 2013 disponibles en la página web de la Superintendencia de Bancos e Instituciones Financieras.

El siguiente código carga el archivo .csv a la variable default, luego asigna cada valor a la región correcta en el objeto gadm, segmenta la variable en cortes, y finalmente asigna las etiquetas apropiadas para que se visualice mejor.

default <- read.csv(file='default.csv', sep = ';', header = TRUE)
gadm$Default <- sapply(gadm$NAME_1, 
                       function(x){default$Default[x == default$Region]})
gadm$Default_Cut <- as.factor(cut(gadm$Default, 
                                  seq(from = 0.02, to = 0.075, by = 0.011)))
levels(gadm$Default_Cut) <- c("<3,1%", "3,1-4,2%",
                              "4,2-5,3%", "5,3-6,4%", "6,4-7,5%")

Paso 5: ¡Graficar!

Ahora graficamos. Lo primero es crear una paleta de colores. El paquete RColorBrewer tiene muchas paletas distintas, de acuerdo a lo que se desee hacer. El paquete distingue tres tipos de paletas de colores: Secuenciales, como una escala de grises, que permiten mostrar escalas de valores; divergentes, como una que parta del rojo, pase por el morado y termine en el azul, resaltando tanto los valores extremos como los intermedios; y las cualitativas, que poseen muchos colores distintos y sirven para mostrar datos que no tienen ninguna estructura.

Nosotros usaremos una paleta secuencial que va del rojo al verde, y lo invertiremos para que corresponda con nuestros valores (rojo para default alto y verde para default bajo).

Colores <- brewer.pal(5, "RdYlGn") # Cinco colores, uno para cada nivel.
Colores <- rev(Colores) # Invierte paleta.

Con esto terminado, podemos graficar. Un único comando genera el gráfico buscado.

spplot(gadm, "Default_Cut", col.regions = Colores, 
       main = 'Porcentaje de Stock en Default por Región (SBIF)')

El gráfico final queda así:

Stock en Default

Créditos (stock) en default según SBIF a Marzo 2013.

¡Y listo! Hemos creado un mapa de Chile con datos propios en R. En el siguiente post mostraré cómo realizar esto usando Shapefiles y los archivos oficiales de la Biblioteca del Congreso Nacional.

Edit: Actualizados los links al archivo de Default y al shapefile con los mapas en R. ¡Gracias @gonzacisternas!