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,

library( mvtnorm )
library( MASS )

n <- 100
sigma <- 0.5

c.r <- data.frame( x = -1:1, y = 1:-1 )
c.v <- data.frame( x = c(-1,1), y = c(-1,1) )

muestra <- function( n, centros, sigma ){
	n.x.centro <- sample( nrow( centros ), n, replace = T )
	tmp <- mvrnorm( n, mu = c(0,0), Sigma = diag( 2 ) * sigma )
	tmp <- tmp + centros[n.x.centro,]
	tmp
}

muestra.r <- data.frame( clase = "red", muestra( n, c.r, sigma ) )
muestra.v <- data.frame( clase = "green", muestra( n, c.v, sigma ) )

mi.muestra <- rbind( muestra.r, muestra.v )

plot( mi.muestra$x, mi.muestra$y,
      col = as.character( mi.muestra$clase ) )

es decir, así:

Los puntos rojos están distribuidos según , una mezcla de tres distribuciones normales esféricas con centros en los puntos (-1,1), (0,0) y (1,-1) y desviación estándar 0.5. Los verdes, según , una distribución similar aunque con centros en (-1,-1) y (1,1):

veros <- function(w, medias, sigma = 5){
	mean( dmvnorm( medias, w, sigma = diag( length( w ) ) / sigma ) )
}

veros.malla <-function(malla, medias, sigma = 5){
	apply( malla, 1, function( x ) veros( x, medias, sigma ) )
}

tmp <- seq( -3, 3, by = 0.03 )
malla <- cbind( x = rep( tmp, each = length(tmp) ),
                y = rep( tmp, length(tmp)) )

bayes.r <- veros.malla( malla, c.r, sigma )
bayes.v <- veros.malla( malla, c.v, sigma )

contour(x=tmp,y=tmp,z=matrix(bayes.r, length(tmp)) )
points(muestra.r$x, muestra.r$y,col="red")

contour(x=tmp,y=tmp,z=matrix(bayes.v, length(tmp)) )
points(muestra.v$x, muestra.v$y,col="green")

Gráficamente,

y

La frontera bayesiana óptima está basada en el siguiente razonamiento, valga la redundancia, bayesiano:

donde:

• es la probabilidad de que la clase/color sea C en el punto x.
• es la probabilidad de observar el valor x.
• es la probabilidad de que x sea un punto de la clase/color C.
• es la probabilidad de la clase C, que suponemos igual en nuestro ejemplo (aunque sería muy desigual en un problema de detección del fraude, por ejemplo, donde el porcentaje de casos fraudulentos es muy bajo).

Bajo las condiciones anteriores,

y, no habiendo penalizaciones asimétricas según la dirección del error, el criterio óptimo de clasificación es asignar la clase R cuando y, de acuerdo con la fórmula anterior, cuando . Es decir, la frontera entre las zonas en que es más probable que una observación proceda de una u otra distribución (es decir, sea de una u otra clase) es

contour(x=tmp,y=tmp,
        z=matrix(bayes.v - bayes.r, length(tmp)),
        levels = 0, labels = "" )
points( mi.muestra$x, mi.muestra$y, col = as.character(mi.muestra$clase) )

que tiene, para este caso concreto, el siguiente aspecto:

Y mañana, ¡a Iguazú!

Escribo desde mi retiro vacacional, en el hemisferio inhabitual, sin wifis y casi de memoria para completar la historia que comencé hace dos semanas en esta bitácora.

Tropecé con el juego que describí en el libro A Mathematician Plays The Stock Market, de John Allen Paulos. Y creo que se equivoca en las probabilidades de los juegos: si en lugar de las que indiqué en mi primera entrada utilizo las suyas, me da la impresión de que el tercer juego es perdedor. ¿Será un bug en el libro? (¿O es que la dislexia me volvió a confundir?)

Paulos cita a el trabajo del físico español Juan Parrondo , en cuya página hay una discusión muy accesible sobre estos juegos paradójicos así como artículos algo más sesudos sobre el tema.

Y más allá de las referencias y las debidas gentilezas con respecto a mis fuentes, aprovecho de pasada para recordar a mis lectores el concepto de supermartingala y el que podría ser considerado su teorema más importante.

Una martingala es un proceso aleatorio tal que . Y una supermartingala, uno en el que .

¿Extraño? Supóngase que es la cantidad acumulada en el juego propuesto. La expresión es el promedio de las posibles cantidades acumuladas tras jugar partidas habida cuenta del resultado de la n-ésima. Por ejemplo, si al cabo de 30 partidas has acumulado 12 euros, puede ser 13 euros con probabilidad 0.49 u 11 con probabilidad 0.51 y, por lo tanto, .

De la discusión anterior se deduce que el primer juego es una supermartingala y, si se jugase con una moneda no sesgada, sería una martingala.

Y el teorema fundamental de las martingalas dice que no hay estrategia capaz de vencer a la banca si esta te ofrece participar en un juego que es una supermartingala (como lo es la ruleta, por ejemplo).

Una posible estrategia para vencer al casino sería jugar aleatoriamente a dos juegos distintos, por ejemplo, la ruleta y el nosequé (no soy aficionado a los divertimentos aleatorios: no sé de otros juegos). El teorema demostraría la imposibilidad también de vencer a la banca usando esta estrategia.

¿Y qué de nuestros juegos y, en particular, el tercero? ¿Cómo puede ser que una estrategia que consiste en alternar entre dos juegos perdedores (dos supermartingalas) nos permita vencer a la banca? ¿Fallan las matemáticas?

Por supuesto que no. Y es porque el segundo juego no es una supermartingala: cuando n no es múltiplo de tres, no se cumple .

El tercero de los juegos, como consecuencia, es uno de esos ejemplos en que se viola ligeramente las condiciones para que se cumpla el teorema para llegar a una conclusión opuesta a él.

Cierro esta discusión con dos ejercicios para mis lectores:

• Probar que con las probabilidades que se indican en el libro de Paulos, el ejemplo no funciona, i.e., que hay un error en su libro. ¡O no!
• Identificar el otro gran contraejemplo del del teorema del tiempo de espera que se cita por doquier y la ciudad rusa a la que se asocia.

Hablamos ya hace un tiempo de las micromuertes. Ahora toca traer a la atención de mis lectores un concepto asociado, el de las microvidas.

Una microvida corresponde a una esperanza de vida de media hora. Malgasta una microvida quien fuma dos cigarros, bebe siete unidades de alcohol (equivalentes a  un litro de cerveza) o vive un día con un sobrepeso de 5 kg.

Microvidas y micromuertes son conceptos análogos, pero no enteramente equivalentes. Ambos nos ayudan a cuantificar pequeños riesgos. Sin embargo, el efecto de las microvidas es acumulativo mientras que el de las micromuertes no: quien haya terminado vivo su sesión de parapente, habrá puesto a cero su contador de micromuertes, pero no así quien haya fumado su segundo cigarro.

En el artículo original de Spiegelhalter aparece también un cálculo aproximado del valor de una microvida y una micromuerte:

El UK National Institute for Health and Clinical Excellence (NICE) sugiere que el National Health Service pague hasta 30.000 libras por un tratamiento que alargue la vida en un año, es decir, unas 17.500 microvidas. Esto equivale a valorarlas en unos £1.70. El UK Department of Transport valora una vida estadística en 1.600.000 libras. Esto significa que están dispuestos a pagar 1,60 libras para evitar una probabilidad de muerte de una en un millón, es decir, una micromuerte.

¿Serán similares las cifras en España? ¿Alguna idea sobre cómo averiguarlo?

Es algo viejo, pero vale la pena traerlo a estas páginas. Se trata de un manifiesto que comienza parejo a aquel otro ahora arrumbado: Un espectro recorre los mercados — el espectro de la falta de liquidez, la congelación del crédito y el fracaso de los modelos financieros.

Habla, sí, principalmente, de finanzas. Pero en gran medida desde la óptica de la modelización y de su responsabilidad en el caos que vivimos ahora. Y, aunque no tiene desperdicio, su colofón de es de universal aplicación y provecho para los modelizadores todos, incluidos los ajenos al mundo de las finanzas. Es una suerte de juramento hipocrático para modelizadores con las siguientes cinco promesas:

• Tendré presente que el mundo no es obra mía y que no satisface mis ecuaciones.
• Aunque utilizaré modelos para estimar valores, no me dejaré influenciar excesivamente por las matemáticas.
• Nunca sacrificaré la realidad en aras de la elegancia sin explicar el motivo.
• Tampoco crearé en los usuarios de los modelos falsas expectativas de precisión; en lugar de eso, haré explícitas las hipótesis y simplificaciones.
• Entiendo que mi trabajo tiene un impacto profundo en la sociedad y la economía que, en gran medida, me es desconocido.

He decidido limpiar mi cartera. Llevo en ella unos cuantos artículos impresos que me acompañan desde hace mucho y que, por un lado, me da pena tirar y, por el otro, no me aportan en el día a día. Voy a reciclar el papel sobre el que los imprimí y, a la vez, dejar en enlace a ellos por si a mí un día (o a alguno de mis lectores otro) me da por volver sobre ellos. Son:

• Don’t be loopy: Re-Sampling and Simulation the SAS Way, del que ya he hablado en otras ocasiones.
• Dealing with Outliers: Maintain your Data’s Integrity, un estupendo resumen sobre el origen, naturaleza y efectos de los outliers. No hay que perderse tampoco algunos de las referencias.
• A practical solution to the pervasive problems of p values, que trata algunos problemas, como este, que plantea el uso de los p-valores y propone una alternativa bayesiana como remedio.
• Is Portfolio Theory Harming Your Portfolio?, un artículo en el que se examinan conceptos que la teoría moderna de los mercados financieros asume acríticamente. Por ejemplo, el autor se rebela frente a la asimilación de riesgo y varianza, cuando el primero debería ser, netamente, la probabilidad de perder dinero. Y examina los límites reales a la diversificación y la presunta eficiencia de los mercados.
• The code of Nature: making authors part with their programs, que repasa algunos fraudes científicos y discute la necesidad que los resultados publicados sean reproducibles y reconstructibles.
• Statistical Modeling: The Two Cultures, en el que Leo Breiman vuelve sobre la lábil frontera entre la estadística y, por llamarlo de alguna manera, la minería de datos. Imperdibles las réplicas y contrarréplicas de Efron y Cox.

Un data scientist es un señor que sabe de varias cosas que no se enseñan ni juntas ni bien por separado en nuestras universidades. Y que, además, se desaprenden rápido en las oficinas y covachuelas donde acabamos ejerciendo. A no ser, claro está, que uno tenga la vocación y la capacidad para nadar contracorriente.

Extraigo de dataists el siguiente gráfico,

que indica cuáles son los tres elementos técnicos —obviando los pertenecientes a otras dimensiones— fundamentales de los que se nutre una carrera como científico de datos.

¿Cómo calificáis vuestras aptitudes en cada una de las tres grandes áreas? ¿Qué podríais hacer por mejorarlas?

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.

Queremos revisar esta comparación bajo otra perspectiva mostrando cómo SAS es el lenguaje perfecto para manejar estructuras de datos complejas y lo fácil que resulta implementar con él algoritmos complejos.

(Nota: con R basta una línea para leer e imprimir un conjunto de datos delimitado: print( read.table( "fichero.txt", sep = ";" )).)

En particular, muestra un pedazo de código para resolver el problema de las ocho reinas. El problema se reduce a encontrar una permutación de los números 1:n tales que

%Macro FirstOf(List);%Scan(&List,1)%Mend;
%Macro RestOf(List);
  %Local lth;
  %Let lth=%Length(%FirstOf(&List));
  %If %Length(&List)>&lth %Then %Left(%Substr(&List,%Eval(1+&lth)));
%Mend;

%Macro OkToAdd(Element,At=,To=,StartAt=);
  %If &To eq %str() or &Element eq %str() %Then 1;
  %Else %If %Sysfunc(Abs(%Eval(%FirstOf(&To)-&Element)))=
    %Sysfunc(Abs(%Eval(&At-&StartAt))) %Then %Do; 0 %Return;%End;
  %Else
    %OkToAdd(&Element,At=&At,To=%RestOf(&To),StartAt=%eval(1+&StartAt));
%Mend;

%Macro qIter(PartialSolution=,List=,Level=,CounterName=);
  %Local item preFix sufFix;
  %If &List eq %str() %Then %Do;%Let &CounterName=%eval(1+&&&CounterName);
    %Put &&&CounterName [&PartialSolution];
  %End;
  %Else %Do;
    %let preFix=;%let item=%FirstOf(&List);%let sufFix=%RestOf(&List);
    %Do %Until (&preFix eq &List);
      %If %OkToAdd(&item,At=&Level,To=&PartialSolution,StartAt=1) %Then
        %qIter(PartialSolution=&PartialSolution &item,
          List=&preFix &sufFix,
          Level=%eval(&Level+1),
          CounterName=&CounterName
        );
        %let preFix=&preFix &item;%let item=%FirstOf(&sufFix);
        %let sufFix=%RestOf(&sufFix);
    %End;
  %End;
%Mend;

%let c=0;
%qIter(PartialSolution=,List=1 2 3 4 5 6 7 8,Level=1,CounterName=c)

Pero me he entretenido en implementar el mismo algoritmo con R y he aquí el resultado:

perm <- function( p, l ){
  foo <- function( x )
    ( a <- length( p ) ) == 0 || all( abs( a:1 ) != abs( l[x] - p ) )
 
  if ( length(l) == 0 )
    cat( p, "\n" )
  else
    invisible( sapply( Filter( foo, 1:length(l) ),
      function( i ) perm( c( p, l[i]), l[-i] ) ) )
}
 
perm(c(),1:8 )

No hay más color que el del resaltador de sintaxis, creo. Y en cuanto a la introducción del artículo, serán mis lectores los que habrán de decidir si tiene más que ver con las nueve que con las ocho reinas.

Informaba el otro día en estas páginas sobre la formación en captura de datos que hubo, dentro del ciclo charlas sobre de periodismo de datos, en MediaLab Prado.

Al final, una de esas imprevisibles circunstancias de última hora me impidió asistir.

Pero, menos mal, ya están disponibles los vídeos para echarles un vistazo este fin de semana.

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!

¿Es posible eso? No, obviamente… ¿o sí? Sigamos leyendo.

En la entrada original proponía tres juegos. El primero, descrito con código así,

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 )

es simple: se tira una moneda y si sale cara, recibes un euro y, si sale cruz, lo pierdes. Aunque la moneda tiene un sesgo de manera que se gana/pierde con probabilidad 0.49/0.51. Como puede observarse, en el largo plazo se tiende a perder dinero si se juega repetidamente.

El segundo juego es parecido al primero pero algo más complejo:

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 )

Se juega con dos monedas:

• si la cantidad ganada hasta la fecha es múltiplo de 3, se juega con una con la que la probabilidade de ganar es de sólo el 9.5 % pero
• en el caso contrario, se juega con otra con la que la probabilidad de ganar es del 74.5 %.

Este es otro juego en el que, jugando repetidamene, también se acaba perdiendo. Para verlo, sólo hay que darse cuenta de que las situaciones en que la cantidad ganada o perdida es múltiplo de 3 representan una especie de barrera probabilística: en ellas casi siempre se pierde. Y las probabilidades de ganar y perder con ambas monedas se han elegido de tal manera que es —un poquito — más probable pasar de tener 3n euros a 3(n-1) euros que a tener 3(n+1) euros.

Así, si observamos el juego sólo cuando la cantidad acumulada es múltiplo de tres (obviando lo que pasa en las jugadas intermedias), veremos que se trata de un juego parecido al anterior.

El tercero de los juegos es más interesante. Es similar a los anteriores, sólo que utiliza una moneda más. En cada jugada, se tira la última moneda y, si sale cara (con probabilidad 0.5), se juega al primero de los juegos y, si sale cruz, al segundo:

juego.fin <- function( x ){
  sample( c( juego.c, juego.s), 1 )[[1]](x)
}

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

hist( res.juego.fin )
fivenum( res.juego.fin )

Lo sorprendente de este juego es que, con él, ¡se gana dinero! Es decir, jugando al azar a uno u otro juego perdedores, puedes acabar jugando a un juego ganador.

¿No es increíble que el bar pueda acabar ganando dinero perdiéndolo en cada venta?

La semana que viene, os cuento el rollo matemático que hay detrás de este singular fenómeno. ¡Prometido!

(Quizás, entre tanto, alguien quiera ir avanzando y enterándose de qué es una supermartingala y de qué va eso del teorema del tiempo de parada óptimo que el ejemplo anterior parece violar).

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).

¡Felicidades a todos los usuarios y entusiastas de R!

Nota: he olvidado mencionar que debo la noticia a Raúl Vaquerizo, que, digámoslo así, me la ha cedido en primicia.