|
Archivo
Archivo para la categoría ‘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,
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
que tiene, para este caso concreto, el siguiente aspecto:
Y mañana, ¡a Iguazú!
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
El código de SAS con el que resuelven este problema es así de estético, expresivo y comprensible:
%Macro FirstOf(List);%Scan(&List,1)%Mend;
%Macro RestOf(List);
%Local lth;
%Let lth=%Length(%FirstOf(&List));
%If %Length(&List)><h %Then %Left(%Substr(&List,%Eval(1+<h)));
%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:
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.
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.
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  independientes, los m puntos del espacion n-dimensional  siguen una distribución uniforme sobre su esfera (su superficie, vale la pena reiterar) unidad.
Para muestrear la bola n-dimensional, hay que muestrear primero la esfera (como en el párrafo anterior) y luego generar m variables aletorias  con la distribución uniforme. La muestra en la esfera unidad será entonces $U_i^{1/n} X_i/\| X_i \|$.
Efectivamente, proporciona la dirección. Y en cuanto a la distancia con respecto al centro hay que tener encuenta que la bola de radio r < 1 contiene sólo un  del volumen de la bola unidad. Como  , entonces  .
En R,
Y si alguien quiere ver las rodajas de una distribución uniforme sobre la esfera 5-dimensional, por ejemplo, puede ejecutar
¿Le sugiere el gráfico (p.e., los rangos de las proyecciones sobre los distintos ejes) algún comentario?
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:
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
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 )
Es una cosa tan maravillosa que no les voy a robar la oportunidad de averiguar por sí mismos en qué consiste. La semana que viene, en la segunda entrega, comentaré el código anterior y explicaré a qué se refiere y, si nadie lo ha dado a conocer antes, dónde reside lo miraculoso del asunto.
Y en la tercera, la última, ofreceré referencias y plantearé un problema que no sé si alguien podrá resolver (y si lo resuelve en la dirección que sospecho, habrá encontrado un sorprendente bug donde menos cabría esperarlo).
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.
Y todo en tres meses. O menos.
Editado en HTML, PDF, etc. Incluso en papel. Con su ISBN y todo. Creative Commons, por supuesto. Con el código disponible en línea en alguna parte.
¿Quimera?
La manera:
- Un coordinador global de la obra, dividida en varios capítulos.
- Un coordinador para cada capítulo, dividido en varias secciones.
- Un redactor para cada sección, consistente en 3-4 hojas que traten un problema concreto.
El coordinador del capítulo es responsable de garantizar cierto grado de coherencia (estilística, etc.) intracapítulo. El coordinador global es responsable de garantizar la coherencia intercapítulo.
Sólo hacen falta 30-40 voluntarios.
Insisto, ¿quimera?
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.
Imaginemos una hipoteca de 100k euros a 25 años referenciada al euribor (2%) y con un diferencial del 1%. Según una calculadora de hipotecas, la cuota es de 474.21 euros mensuales. Y, efectivamente, 25 * 12 * 474.21 = 142263 excede con mucho la cuantía inicial.
Pero —cosa que nunca nos enseñaron a los matemáticos— dinero de hoy no puede sumarse con dinero de mañana. Son peras y manzanas. Cinco euros mañana no son cinco euros en la mano. Los cinco euros de mañana —por muchos motivos: por la inflación, porque quien me los tiene que entregar a lo peor se olvida, porque los necesito ya, etc.— valen menos que los cinco que tengo. Digamos que el valor de cinco euros mañana son, por ejemplo, 4.98 euros hoy porque me veo obligado a aplicarles un descuento.
Los contables, cuando tienen que sumar cantidades en fechas distintas, utilizan el descuento de todas ellas a día de hoy —conocido como valor presente neto— para agregarlas. El descuento que aplican viene determinado por los tipos de interés: 5 euros en 3 meses valen  hoy, donde  es el tipo de interés mensual.
Para calcular el valor de la hipoteca, se descuentan los pagos por la tasa mensual, que en nuestro caso es de 3 / 12 = 0.25 % y se agregan de la siguiente forma:
Aunque la expresión anterior admite una forma cerrada, nos es más cómodo a los perezosos calcularla en R así
> sum( ( 1 + 0.03 / 12) ^(-(1:(25 * 12)) ) ) * 474.21
[1] 99999.72
Y, efectivamente, se obtienen los 100k euros de partida.
Pero, entonces, ¿dónde y cómo gana dinero el banco?
Pues resulta que no todos los agentes económicos descuentan con los mismos tipos de interés. El Estado Español tiene unos tipos de descuento distintos que los del Alemán o que la Comunidad Valenciana. Y mi banco otros distintos que los míos: le cuesta menos endeudarse que a mí.
Supongamos que el banco que concede la anterior hipoteca es capaz de financiarla pagando el euribor, es decir, el 2%. ¿Cuál sería el valor presente neto de los flujos de dinero que supone el pago de las cuotas? Es
que en R queda
> sum( ( 1 + 0.02 / 12) ^(-(1:(25 * 12)) ) ) * 474.21
[1] 111880.4
Et voilá, el banco hace casi 12k euros con conceder dicha la hipoteca.
|
