R

Una de las preguntas formuladas dentro del foro desde el que seguimos la lectura del libro The Elements of Statistsical Learning se refiere a cómo construir la frontera bayesiana óptima en ciertos problemas de clasificación.

Voy a plantear aquí una discusión así como código en R para representarla (en casos simples y bidimensionales).

Supongamos que hay que crear un clasificador que distinga entre puntos rojos y verdes con la siguiente pinta,

No sé si habéis visto la película argentina Nueve reinas. Trata de unos timadores que engatusan a incautos para sacarles la platica.

Pero no voy a hablar de esas nueve reinas sino de las ocho de Solve Eight Queens Puzzle With SAS Macro. De su introducción extraigo y traduzco:

The Little SAS Book contiene un excelente ejemplo para ilustrar las diferencias entre SAS como lenguaje de programación y C++ mostrando lo complicado que puede resultar procesar conjuntos de datos con un lenguaje de propósito general. Son 28 líneas de código C++ y 5 de SAS para leer un fichero delimitado e imprimirlo por pantalla. Es un ejemplo perfecto de cómo SAS es un lenguaje de cuarta generación con un alto nivel de abstracción y expresividad.

Tal como prometí hace ahora una semana, voy a añadir las palabras que faltaban en aquella entrada. Pero primero, imaginad un bar en el que se venden cafés y cervezas. El coste de servir un café es de 1.10 euros pero se vende por 1. El coste de servir una cerveza es 1.30 euros pero se vende por 1.10. Entran los clientes y piden o café o cerveza. ¡Y resulta que a fin de mes el bar hace dinero!

Más evidencias sobre la emergencia de R: ha entrado en el top 20 de lenguajes de programación elaborado por Tiobe por primera vez en enero de 2012:

La lista, según avisa el mismo Tiobe, no es científica: se basa en un índice de popularidad elaborado a partir de información de ofertas laborales, buscadores de internet, etc.

Nótese también que es el primero de los lenguajes de programación que no es de propósito general sino de dominio (la estadística). Nótese también su posición relativa con respecto a Matlab (23) y SAS (32).

Debo esta entrada a la diligencia de Juanjo Gibaja, que se tomó la molestia de ubicar los teoremas relevantes en el libro Simulation and the Monte Carlo Method de Rubinstein y Kroese.

Esencialmente, como la distribución normal multivariante (con matriz de covarianzas I) es simétrica, entonces, dadas $latex X_1,\dots, X_m \sim N( 0, I_n )$ independientes, los m puntos del espacion n-dimensional $latex X_i/| X_i |$ siguen una distribución uniforme sobre su esfera (su superficie, vale la pena reiterar) unidad.

Ante las preguntas de alguno de mis lectores, voy a proporcionar una explicación acerca de la misteriosa L. Bueno, voy más bien a dejar que la deduzcan ellos mismos a partir de la siguiente serie de bloques de código:

a <- rep( 0, 10 )
typeof( a )
object.size( a )

b <- rep( 0L, 10 )
typeof( b )
object.size( b )

##############

a <- 1:10
typeof( a )
object.size( a )

a[1] <- 10
typeof( a )
object.size( a )

a <- 1:10
a[1] <- 10L
typeof( a )
object.size( a )

##############

a <- 1:10
tracemem( a )
a[1] <- 2

a <- 1:10
tracemem( a )
a[1] <- 2L

##############

system.time( replicate( 1e5, { a <- (1:100); a[1] <- 12  } ) )
system.time( replicate( 1e5, { a <- (1:100); a[1] <- 12L } ) )

Lectores míos, no seáis perezosos y haced, cuando menos, ?tracemem en vuestra consola. Una vez leída la página de ayuda, ¿se os ocurre algún truco para ahorrar mucha memoria cuando trabajáis con objetos (p.e., matrices) grandes de enteros?

Hoy voy a hacer mención a una cosa prodigiosa. Pero sin palabras. Voy a regalar a mis lectores tres pedazos de código que son este

jugar <- function( n, make.step ){
  tmp <- rep( 0L, n)
  for( i in 2:n )
    tmp[i] <- make.step( tmp[i-1] )
  tmp
}

juego.s <- function( x, prob.perder = 0.51 ){
  x + ifelse( runif(1) < prob.perder, -1L, 1L )
}

res.juego.s <- replicate( 1000, jugar( 1000, juego.s )[1000] )
hist( res.juego.s )
fivenum( res.juego.s )

este

juego.c <- function( x ){
  prob.perder <- ifelse( x %% 3 == 0, 0.905, 0.255 )
  juego.s( x, prob.perder )
}

res.juego.c <- replicate( 1000, jugar( 1000, juego.c )[1000] )

hist( res.juego.c )
fivenum( res.juego.c )

y este otro

No sé si es quimera o no. Se me ocurrió el otro día. Dejo mi idea aquí escrita por ver por dónde respira la comunidad.

Se trata, sí, de un libro extenso sobre R. Que cubra el 90-95% de los métodos estadísticos que utilizan los usuarios —en sentido amplio— de la estadística: médicos, sociólogos, etc. Con R. Con la teoría justa pero, eso sí, con referencias a fuentes fiables: se supone que sus lectores saben ya algo de estadística, pero tal vez no cómo afrontar su problema con R. Una especie de recetario bien hecho. Un vademécum.

Mediante la presente, notifico a los interesados en la lectura de “The Elements of Statistical Learning” que esta semana tenemos que dar cuenta de los capítulos 1 (que es una introducción muy ligera) y 2 (donde comienza el tomate realmente).

Esta noche Juanjo Gibaja y yo estudiaremos la mecánica de lectura en común.

Los interesados pueden escribirme a cgb@datanalytics.com para, de momento, crear una lista de correo.

Parece mentira, pero hay gente que lo calcula fatal. Hace tiempo, un antiguo colega mío, matemático él, había propuesto el ejercicio a sus alumnos y estimó, me contó, que el banco recibía, aproximadamente, el doble de lo que prestaba. La operación que había realizado era muy sencilla: calcular el saldo vivo inicial con la suma de todas las cuotas mensuales. Pero la operación es incorrecta. Veamos por qué. Y obtengamos, de paso, alguna estimación más ajustada.

Un gráfico de pares de variables —que no he sabido traducir mejor desde el original inglés pairplot— es algo como lo siguiente:

Es posible ahora construir gráficos de pares más sofisticados e informativos usando el paquete GGally de R. Usando el código (extraído de SAS and R)

library(GGally)

ds <- read.csv("http://www.math.smith.edu/r/data/help.csv")
ds$sex <- as.factor( ifelse(ds$female==1, "female", "male") )
ds$housing <- as.factor( ifelse(ds$homeless==1, "homeless", "housed") )
smallds <- subset(ds, select=c("housing", "sex", "i1", "cesd"))

ggpairs(smallds,
        diag=list(continuous="density", discrete="bar"),
        axisLabels="show")

se obtiene la siguiente versión mejorada:

Puede que algunos de mis lectores sepan que el lucero del alba es el nombre con que se conoce al planeta Venus cuando es visible en el cielo al amanecer.

En contextos menos poéticos se conoce por tal nombre a esto:

Es decir, una determinada configuración de los precios de apertura y cierre de tres días de cotización (bursátil, por ejemplo) de forma que:

• El primer día hay una bajada
• El tercer día hay una subida
• Los precios de apertura y cierre del segundo día son inferiores a los del cierre del primero y apertura del segundo.

Se ve que eso es cosa güena. De El Economista extraigo el siguiente párrafo atribuido a un tal Joan Cabrero:

Propone Juan José Gibaja como propósito intelectual para el año nuevo el leer The Elements of Statistical Learning —libro que puede descargarse gratuita y legalmente del enlace anterior— de tapa a tapa, en grupo y a razón de capítulo por semana.

La idea es hacerlo en común, enlazando el contenido del libro con código —sea disponible o de nuevo cuño cuando la situación lo requiera— y haciendo públicos las ideas que resulten de esta lectura en una red de bitácoras (a la que esta pertenecería).