|
Archivo
Archivo para la categoría ‘r’
Acaba de terminar la primera reunión del grupo de usuarios de R de Madrid. No hemos disfrutado de la más primaveral de las tardes. Y la ubicación era un tanto excéntrica. Pero hemos tenido tres charlas muy interesantes (y luego, la mía, claro), nueve asistentes (¡espero haber contado bien!) y, sobre todo, unos intercambio de ideas sumamente provechosos.
Los enlaces a las presentaciones estarán pronto disponibles en la página del grupo. Pero como delanto:
La primera, de Rafael Garcia Leiva, de entropycs, sobre el uso de paquetes para el análisis financiero de R aplicados al trading algorítimico, es decir, a ganar dinero por las noches, mientras uno duerme, en los mercados.
Carlos Ortega ha aplicado sus conocimientos de la ingeniería del software al análisis cuantitativo de la calidad del desarrollo de R como proyecto de software a través del historial de bugs a través de las distintas versiones desde 1998.
Finalmente, Pedro Concejero nos ha hecho una demostración en vivo del uso de los paquetes de la familia bigX (bigmemory y relacionados) para el tratamiento y análisis de ficheros de, en este caso concreto, 10GB. También nos ha anunciado que nos contará algo que tiene que ver con colores… ¡y hasta ahí puedo leer!
(Y no quiero hablar de la mía, que era la más inacabada e inconexa. Como imperfecta la he definido horas antes de darla. De hecho, el código que he enseñado era del día de antes y funcionaba a trancas y barrancas. Pero, al fin y al cabo, hablaba de a dónde quería hacerlo llegar. Y esa meta sí que será, si Dios quiere, pública.)
Obviamente, habrá una próxima. Habrá que ver si MediaLab Prado o Tabacalera nos acogen. Porque creo que una ubicación más céntrica sería preferible, ¿no?
Hoy ando demasiado ocupado para escribir. Y como es posible que alguno de mis lectores no lo esté tanto como para no leer, le dejo un artículo de Tukey (abajo del todo en el enlace anterior) para que conozca al personaje, si no ha tenido el gusto previo, disfrute en cualquer caso y, en todos, sepa de dónde vienen los rootograms que implementa el paquete latticeExtra de R.
El miércoles 21 de marzo de 2012, en el aula N-130 del edificio de primer curso (también conocido como Prefabricado) de la facultad de CC. Económicas de la UCM (Somosaguas) tendrá lugar la primera reunión del grupo de usuarios de R de Madrid.
Contamos con tres charlas muy interesantes y una mía. Esta última trata de una función que aún no existe sino en forma de bosquejo en mi cabeza. Espero que esté presentable el miércoles.
De todos modos, espero que el plato fuerte del programa sea el social.
Y quiero animar a los entusiastas de R de Madrid, que me consta que somos muchos, a pasarse por Somosaguas (ya lo sé, una ubicación un tanto excéntrica).
¿Os animáis?
Esta mañana casi me da esa tontería de sentirme orgulloso de ser de donde soy, Zaragoza. Al fin y al cabo, podría haber sido de cualquier otro lugar. Pero es que Zaragoza tiene uno de los portales de datos públicos municipales más avanzados. En eso es una ciudad pionera.
(Se lo hemos de agradecer a nuestro alcalde, Belloch, que, dicen las malas lenguas, además de socialista y barbudo, es linuxero).
Entre los datos disponibles, los hay de tráfico en tiempo real. En particular, existe una serie de tramos de calle y un fichero que se actualiza cada pocos segundos que indica el estado del tráfico en ellos.
Y he pensado que tal vez podría hacer una virguería con R.
Así que he escrito lo siguiente:
library( rjson )
# tmp <- readLines( "http://www.zaragoza.es/trafico/estado/tramos23030.json" )
tmp <- readLines( "http://www.zaragoza.es/trafico/estado/tramoswgs84.json" )
tmp <- fromJSON( tmp )[[1]]
status <- fromJSON( readLines( "http://www.zaragoza.es/trafico/estado/estado.json" ) )
status.time <- status$timestamp
status <- strsplit( status$estados, "" )[[1]]
# length( kkk )
tmp <- lapply( tmp, function( x ) {
id <- x$id
name <- x$name
status <- status[id]
lat <- sapply( x$points, function( y ) y$lat )
lon <- sapply( x$points, function( y ) y$lon )
data.frame( id = id, name = name, status = status, lat = lat, lon = lon )
})
tmp <- do.call( rbind, tmp )
# tmp <- merge( tmp, status )
plot( range( tmp$lon ), - range( -tmp$lat ),
xaxt = "n", yaxt = "n", type = "n",
main = paste(
"Estado del tráfico en Zaragoza",
strptime( gsub( "-|Z", " ", status.time),
format = "%Y%m%d %H%M%S" ), sep = "\n" ),
xlab = "", ylab = "" )
foo <- function( x, y, status ){
colores <- c( "black", "red", "yellow", "green", "lightgray" )
color <- colores[ match( status, c( "b", "r", "y", "g" ), nomatch = 5 ) ]
lines( x,y, col = color, lwd = ifelse( status == "-", 1, 2 ) )
}
by( tmp, tmp$id, function( x ) foo( x$lon, x$lat, status = x$status ) )
Que da como resultado (a la hora en la que lo he ejecutado, cuando los zaragozanos están ya casi todos en su casa)
Pero me ha sabido a poco y he querido hacerlo todavía más a lo maño. Así que he añadido
library( OpenStreetMap )
map <- openmap( c( max(tmp$lat), min(tmp$lon) ), c( min( tmp$lat ), max(tmp$lon) ), type = "osm")
plot(map,raster=TRUE)
tmp.mercator <- data.frame( projectMercator( tmp$lat, tmp$lon ) )
tmp.mercator$status <- tmp$status
foo <- function( x, y, status ){
colores <- c( "black", "red", "yellow", "green", "lightgray" )
color <- colores[ match( status, c( "b", "r", "y", "g" ), nomatch = 5 ) ]
lines( x,y, col = color, lwd = ifelse( status == "-", 1, 2 ) )
}
by( tmp.mercator, tmp$id, function( x ) foo( x$x, x$y, status = x$status ) )
Y he obtenido
Hay algunas cosas que me gustaría poder añadir, minucias, pero que estoy demasiado ocupado para investigar y que me gustaría dejar de tarea a mis lectores:
- ¿Cómo poner un título al segundo gráfico?
- ¿Cómo difuminar la imagen de fondo para que resalten más los tramos de tráfico sobre el excesivo detalle del mapa subyacente?
Que conste que soy un partidario de los adjetivos. Supongo que por sentimentalismo. Me caen simpáticos excepto
- cuando se abusa de ellos y se dice, por ejemplo, analítica en lugar de análisis o normativa en lugar de norma o
- los usan estadísticos en horario laboral.
Y si trabajan en el INE, aún más: se les paga por estadísticos, no por guionistas de opereta.
Viene esto al siguiente párrafo (con mi subrayado):
En el año 2010 murieron a manos de sus parejas o exparejas 73 mujeres, que ejercieron sobre ellas violencia de género. Esta cifra, tras el esperanzador descenso experimentado en el año 2009 (55 mujeres muertas) alcanza el nivel de años anteriores, 76 mujeres en 2008 y 71 mujeres en 2007.
Que puede consultarse en la página dedicada a Víctimas mortales por violencia de género del INE. Por si cambiase el enlace, se accede a dicha página navegando a la sección Delito y Violencia de la publicación gratuita Hombres y mujeres en España. La serie temporal que lo acompaña es:
Y los datos correspondientes, dat <- c( 63, 50, 54, 71, 72, 57, 68, 71, 76, 55, 73 ).
Hay que notar que algunos de estos números corresponden a hombres asesinados por sus parejas (mujeres), tal y como indica el pie del gráfico en el INE, aunque el comentario que lo acompaña y que reproduzco parece ignorar (tal vez excusablemente) este pequeño matiz.
Hay que advertir también que sólo el 0.029% de las mueres fallecidas en el 2009, 55, lo fueron por esta causa, mientras que por, por ejemplo, complicaciones de la atención médica y quirúrgica murieron 241; por gripe, 141; por tuberculosis, 143; por SIDA, 226; por caídas accidentales, 892; por suicidio, 763 y por muerte súbita infantil, 33. Lo cual lleva a uno a preguntarse por las causas del desigual impacto de cierto tipo de noticias en los medios. Pero esa es otra cuestión.
Habida cuenta de lo infrecuente del fenómeno, es decir, el asesinato de una mujer por parte de su pareja, y de la posible independencia entre este tipo de sucesos, cabe pensar que el número anual de casos sigue una ley de Poisson. De ser así, su parámetro (estimado por máxima verosimilitud) sería mean(dat), es decir, 64.54. Dado que el tamaño de la muestra es tan pequeño, no se me ocurre ningún procedimiento para evaluar la bondad del ajuste. No obstante,
Es decir, bajo el modelo propuesto, la mediana del número máximo y mínimo anual de este tipo de asesinatos coincidiría con los observados. No sólo es esto evidencia —heterodoxa y discutible— a favor del modelo sino, también, de cómo el adjetivo esperanzador que usa el INE no tiene tanto que ver con la esperanza (otro de los nombres de la media) sino con la variación.
La variación en las cifras entraría, por tanto, dentro de lo normal —¿de lo Poisson?— y ni habría que entregar medallas en el 2009 ni quitárselas en el 2010 a los probos y sacrificados funcionarios encargados de luchar contra esta tan nimia como antiestética lacra social.
A las entradas que he hecho sobre Julia estos últimos días, quiero añadir esta en la que publico mi primer programa en dicho lenguaje.
Me ha dado por reimplementar el programa para realizar un muestreo de Gibbs que aparece en Gibbs sampler in various languages.
Lo primero ha sido instalar Julia, para lo que basta con seguir las instrucciones que aparecen en su página de github. Y aviso: tarda bastante en descargar y compilar todas sus dependencias.
El código de mi programa, para que quede constancia pública de que vengo a ser un pionero en España usándolo, es:
function gibbs(n, thin)
x = 0.0
y = 0.0
for i=1:n
for j=1:thin
x = randg( 3 ) / ( y * y + 4 )
y = randn() / sqrt( 2*x + 2 ) + 1 / (x + 1)
end
println( i, " ", x, " ", y )
end
end
gibbs( 50000, 1000 )
El programa en cuestión, tal cual viene en la página que cito arriba, en R, ha tardado en ejecutarse
$ time Rscript mcmc00.R > data.tab
real 8m31.259s
user 8m29.668s
sys 0m0.948s
Y en Julia,
$ time julia mcmc00.j > data.tab.j
real 0m9.268s
user 0m9.113s s
ys 0m0.144s
Es decir, del orden de 60 veces menos. Es de reseñar que este tipo de algoritmos son de los menos indicados para ser ejecutados con R: son imposibles de vectorizar. No en vano, MCMC (tal como aparece en los nombres de los programas) significa Markov Chain Monte Carlo y en una cadena de Markov, cada valor depende del estado anterior. Hay que morir a hierro y construir un bucle.
El resultado de ambas simulaciones puede apreciarse en el siguiente gráfico:
Es cierto que en este ejemplo concreto, los números aleatorios podrían precalcularse en una única llamada a rnorm o rgamma, pero no he observado una ganancia sustancial de tiempos, que se pierden, fundamentalmene, en el código interpretado dentro del bucle.
No sé si conocéis Julia, un lenguaje de programación orientado al cálculo científico. Os dejaré echarle un vistazo a su página.
¿Ya?
Bueno, pues estoy un poco enfadado con ellos. Me pasa un poco como a los catalanes que se quejaban de que en las fotos de ABC siempre sacaban a Jordi Pujol (todavía más) feo (de lo que por sí era): en las comparaciones no le hacen excesiva justicia a R. Me he tomado la molestia de reescribir el código para una de las comparaciones que realizan, pi_sum, utilizando código vectorizado.
El original es
library(R.utils)
assert = function(bool) {
if (!bool) stop('Assertion failed')
}
timeit = function(name, f, ..., times=5) {
tmin = Inf
for (t in 1:times) {
t = system.time(f(...))["elapsed"]
if (t < tmin) tmin = t
}
cat(sprintf("r,%s,%.8f\n", name, tmin*1000))
}
pisum = function() {
t = 0.0
for (j in 1:500) {
t = 0.0
for (k in 1:10000) {
t = t + 1.0/(k*k)
}
}
return(t)
}
assert(abs(pisum()-1.644834071848065) < 1e-12);
timeit("pi_sum", pisum, times=1)
y tarda 10760 milisegundos en mi máquina.
El código alternativo,
pisum.cjgb <- function() {
for (j in 1:500)
t <- sum( (1:10000)^(-2) )
return(t)
}
assert(abs(pisum.cjgb()-1.644834071848065) < 1e-12);
timeit("pi_sum_cjgb", pisum.cjgb, times=1)
tarda 510 milisegundos, veinte veces menos.
En cualquier caso, me viene sorprendiendo mucho la velocidad de JavaScript. ¿Recordáis mi entrada sobre las ocho reinas? La versión en JavaScript,
<body onload="javascript:cnt=0;
function backTrack(trial,next){
if (trial.length==0){return true;}
else {
for (var i in trial){
if (Math.abs(trial.length-i)==Math.abs(next-trial[i])){
return false;
}
}
return true;
}
}
function perm(p,l){
if (l.length==0){
cnt++;document.write(cnt+'::'+p+'<br/>')
}
else {
for (var i in l){
if (backTrack(p,l[i])){
perm( p.concat(l[i]),
l.slice(0,l.indexOf(l[i])).concat(l.slice(l.indexOf(l[i])+1,l.length))
);
}
}
}
}
perm([],[1,2,3,4,5,6,7,8,9,10,11])">
</body>
(cópiese el código anterior en un fichero loquesea.html y ábrase en el navegador) es notablemente más rápido que el equivalente en R.
¡Supongo que hay mucha más gente optimizando JavaScript —que, al fin y al cabo, es el instrumento que usa Google para entrar a tu cocina— que R!
… cuando dice que hay que ver qué pasa y analizar las estadísticas. En lo demás, no lo sé (ni lo pienso decir aquí). Pero traigo el asunto a colación porque hace un par de días hablé, un tanto exteporáneamente, sobre desempleo y subsidios. Y uno de mis lectores hizo un comentario del que extraigo
No me gusta, no me gusta que se insinúe siempre que “España está llena de listos, que agotan el paro porque les sale mejor que trabajar”.
Ahora, Juan Rosell, presidente de la CEOE, nos dice:
[...] hay que ver qué pasa y analizar las estadísticas. Y estas dicen que como aquí el subsidio dura hasta 24 meses, la gente encuentra trabajo milagrosamente cuando falta un mes o dos para agotar el subsidio. Eso quiere decir que no está funcionando del todo bien.
Obviamente, en estas páginas no podemos estar más de acuerdo con el Sr. Rosell en cuanto a la necesidad de echarle un buen vistazo a los números. Y como él no nos dice más y yo ando con mucho trabajo estos días, aprovecho y pregunto:
- ¿Sabrá alguno de mis lectores a qué estadísticas se refire el Sr. Rosell? ¿Están disponibles?
- ¿Puede probarse estadísticamente su afirmación? ¿Nos puede regalar al resto de los lectores la demostración?
La serendipia me llevó a toparme con el RUGBCN, es decir, el grupo de usuarios de R de Barcelona. Me puse en contacto con ellos y Lluis Ramon ha tenido la gentileza de ofrecerse a responder una serie de preguntas mías que espero que, por un lado, animen a los usuarios de R de BCN a acercarse a las reuniones y, por otro, sirvan de estímulo para la creación de grupos de usuarios similares en otros lugares.
He aquí la entrevista:
¿De quién partió la iniciativa de crear un grupo de usuarios en BCN? ¿Cómo se dio a conocer?
La iniciativa salió de Tim que lo propuso en el Local R User Group Directory de Revolutions. Allí hay un apartado de groups wanted y uno de ellos era Barcelona. En el apartado de los comentarios estaba su mensaje. Cuando yo lo vi llevaba casi un año publicado. Sin dudarlo le envié un correo preguntando si continuaba interesado y, por suerte, lo estaba.
¿Cómo, dónde y cuándo os soléis reunir? ¿Cómo es una reunión típica?
Generalmente nos reunimos los jueves de 19h a 21h cada mes o mes y medio. Al principio en un bar de Barcelona que tiene una sala con proyector y nos la facilitaba. Últimamente nos hemos reunido en una sala muy acondicionada en el PRBB que nos la reserva Maik. Todavía no hay reuniones típicas ya que estamos probando cual es el formato que más nos gusta. La primera fue de presentación y enseñamos 3 IDE’s que usamos: Eclipse, Emacs y RStudio. En la segunda se hicieron dos presentaciones del tipo lección magistral y en la última formamos 3 grupos pequeños, dónde se trataron los temas de forma muy interactiva en la primera hora y en la segunda visualizamos juntos el webinar de Hadley presentando las novedades de ggplot2. Lo único que hay típico es que luego nos vamos unos cuantos a tomar una cervecita y a comer algo.
¿Cuál es el perfil de los asistentes a las reuniones? ¿Tenéis intereses comunes (bioestadística, etc.)? ¿O procedéis de entornos diversos?
El grupo es muy heterogéneo. Desde usuarios que saben mucho al que viene por curiosidad. La mayoría somos locales, pero también hay gente de USA, Alemania, Italia, etc que están trabajando o haciendo doctorados en Barcelona. Aunque mayoritariamente la gente proviene del mundo académico-universitario hay miembros de la administración pública, del sector privado como de alguna start up. Un interés común claro es el R jeje. Claramente abunda la bioestadística pero hay otros campos en el grupo como demografía, tráfico, simuladores, energía, consultoría, diseño de fármacos,..
¿Cuál es la cosa más curiosa, útil o provechosa que has aprendido en alguna de las reuniones?
Me gustó mucho la presentación de Aleix sobre Reference Classes. Sabía que existían pero nunca me había puesto a mirarlo a fondo, tampoco había mucha documentación porque son relativamente nuevas. Como él las usa día a día, fue una introducción mucho más amena que buscando por Internet. También me gustó ver como funcionaba Emacs porque había leído que mucha gente lo usaba.
¿Qué consejos podéis dar a quienes proyecten organizar grupos similares en otros lugares?
Que no se lo piensen mucho y que lo hagan. Que tampoco necesita mucha planificación ni organización, te reúnes con gente que les gusta lo mismo que a ti y habláis del tema. La primera vez en un lugar tranquilo, en un bar y explicáis que es lo que hacéis, que packages preferís de R o que IDE usáis, etc.
Finalmente, ¿nos veremos en las IV Jornadas de R?
Me gustaría mucho, aunque no sé si será compatible con mi trabajo. Lo que estoy seguro es que alguien de RUGBCN sí que irá y nos mantendrá informados de cómo han ido y qué se ha presentado.
A la pregunta, tal vez con una formulación mejorable de un usuario de la lista de R, sobre cómo representar una distribución normal bivariada con correlación 0.5 en 3D di ayer esta solución:
library( mvtnorm )
x <- y <- -20:20 / 10
z <- matrix( 0, length( x ), length( y ) )
m <- c(0,0)
sigma <- matrix( c(1, 0.5, 0.5, 1 ), 2 )
for( i in 1: length( x ) )
for( j in 1:length( y ) )
z[i,j] <- dmvnorm( c( x[i], y[j] ), c(0,0), sigma )
persp( x, y, z )
No obstante, la solución alternativa de Carlos Ortega es toda una virguería que merece ser reproducida en estas páginas:
library(fMultivar)
library(rgl)
x = (-40:40)/10
X = grid2d(x)
z = dnorm2d(X$x, X$y, rho = 0.5)
Z = list(x = x, y = x, z = matrix(z, ncol = length(x)))
open3d()
bg3d("white")
material3d(col="black")
persp3d(Z$x, Z$y, Z$z, aspect=c(1, 1, 0.5), col = "lightblue", xlab = "X",
ylab = "Y", zlab = "Z")
play3d(spin3d(axis=c(0,0,1), rpm=5), duration=20)
¿Os gusta?
|
