España, ¿radial? (II)

Una de las principales objeciones que se le pueden hacer a mi entrada de ayer es que puede estar confundiendo la causa con efecto: puede que parte de la radialidad de la red que obtuve tenga que ver con el tamaño desproporcionado de Madrid que, a su vez, podría haber sido causado por la radialidad de la red tradicional de las comunicaciones españolas. Así que enviemos una partida de pescado en malas condiciones a Mercamadrid, convidemos a toda la provincia, veámosla fenecer víctima de contumaces diarreas y rehagamos la simulación suponiendo que ...

26 de abril de 2012 · Carlos J. Gil Bellosta

España, ¿radial? (I)

Me propuse hace un tiempo combinar lo que aprendí creando rutas callejeras por Zaragoza con una entrada que escribí sobre la estructura radial de las vías de transporte de España. El problema que me planteo es si tiene sentido que la red de carreteras Española tenga estructura radial habida cuenta de la geometría peninsular bajo ciertas hipótesis, siempre discutibles y mejorables, de partida. Así que, en primer lugar, cargué los paquetes de R necesarios, un fichero que creé que contenía las capitales de provincia, su latitud, su longitud y la población de las respectivas provincias y fabriqué una red de carreteras muy ineficiente que unía todos los nodos entre sí: ...

25 de abril de 2012 · Carlos J. Gil Bellosta

Segunda reunión de usuarios de R de Madrid: recordatorio

Aprovecho para recordar a los usuarios de R de Madrid que el jueves 26 de abril, a las siete de la tarde, tendrá lugar la segunda reunión del grupo de usuarios de R de Madrid en la sala Metrópolis de La Tabacalera (glorieta de Embajadores). El programa, como siempre, puede consultarse en la página del grupo.

24 de abril de 2012 · Carlos J. Gil Bellosta

Variables instrumentales con R

Los economistas usan unas cosas a las que llaman variables instrumentales con las que uno apenas se tropieza fuera de contextos econométricos. El problema se plantea en el contexto de la regresión $$y_i = \beta x_i + \varepsilon_i,$$ cuando existe correlación entre X y $\varepsilon$. En tales casos, el estimador por mínimos cuadrados es $$\hat{\beta} =\frac{x’y}{x’x}=\frac{x’(x\beta+\varepsilon)}{x’x}=\beta+\frac{x’\varepsilon}{x’x}$$ y debido a la correlación entre X y $\varepsilon$, está sesgado. La solución que se plantea en ocasiones es el de usar variables instrumentales, es decir, variables correlacionadas con X pero no con $\varepsilon$. La siguiente simulación en R ilustra el problema: ...

19 de abril de 2012 · Carlos J. Gil Bellosta

Rutas por Zaragoza con R

Óscar Perpiñán me puso el otro día al tanto del paquete osmar de R, que proporciona la infraestructura para acceder a datos de OpenStreetMap a través de diferentes fuentes, trabajar con ellos con R de una manera unificada y aprovechando la infraestructura que proporcionan otros paquetes como, por ejemplo, sp e igraph. Hoy voy a ilustrar el uso de este paquete adaptando un ejemplo de sus autores para encontrar la ruta óptima entre dos puntos de Zaragoza, la mercería Bell y el colegio La Salle Montemolín, ambos lugares muy vinculados a mi prehistoria. Comenzaré cargando los paquetes necesarios y los datos de OpenStreetMap correspondientes a Zaragoza: ...

16 de abril de 2012 · Carlos J. Gil Bellosta

Corrección por exposición del modelo logístico

He tropezado con una extensión curiosa y que no conocía del modelo logístico que lo emparenta un tanto con los modelos de supervivencia. Es un problema que aparece en los modelos de los actuarios, por ejemplo, y en la supervivencia de nidos (sí, nidos de bichos alados), parece. Es el siguiente: supongamos que unos sujetos están expuestos a un cierto suceso cuya probabilidad, $p_i$, depende del sujeto a través del esquema habitual de la regresión logística (es decir, depende de algunas variables como el sexo, etc., a través de una fórmula lineal cuyos coeficientes interesa estimar). ...

11 de abril de 2012 · Carlos J. Gil Bellosta

Un intérprete alternativo de R

Java es un lenguaje de programación que puede ejecutarse sobre muchas máquinas virtuales distintas: la de Sun, la de IBM, etc. Algo parecido pasa con SAS, que puede ejecutarse sobre el intérprete de SAS Institute o sobre el de WPS. El código escrito en R puede ejecutarse, en principio, en dos plataformas distintas: La creada por el R Development Core Team y que todos, más o menos, conocemos. La desarrollada por Tibco (y, previamente, por Insightful) para S-Plus, el dialecto propietario de R (o S). ¿Son esas todas las opciones? Sí, por el momento. ...

10 de abril de 2012 · Carlos J. Gil Bellosta

De D'Hondt a Banzhaf

Hablé el otro día con Emilio Torres y comentamos de pasada la situación política en Asturias, donde vive, después de las últimas elecciones. El escaño obtenido por UPyD otorgaba a tal partido un poder en exceso del tamaño de su representación porque era clave para formar el futuro gobierno del principado. Pero, ¿cuánto poder realmente supone ese escaño en esas condiciones? ¿Puede cuantificarse? Porque se habla mucho en periodo electoral de la ley D’Hondt pero, una vez asignados los escaños, cambia el juego. ...

4 de abril de 2012 · Carlos J. Gil Bellosta

Otra de huelgas

Hoy, por motivos evidentes, e igual que en septiembre de 2010, voy a hablar de huelgas. De la misma fuente que entonces he descargado este fichero. Y he ejecutado library(pxR) library(reshape) library(ggplot2) dat <- read.px("pcaxis-623612450.px") dat <- as.data.frame(dat) dat.mes <- cast(dat, Periodo ~ series) colnames(dat.mes) <- c("mes", "n.huelgas", "n.trabajadores", "n.jornadas") p <- ggplot(data = dat.mes) + geom_line(aes(x = mes, y = n.huelgas, group = rep(1, nrow(dat)))) p ggsave("huelgas_por_mes.png") dat.anno <- dat tmp <- strsplit(as.character(dat.anno$Periodo), "M") dat.anno$Periodo <- sapply(tmp, function(x) x[1]) dat.anno <- cast(dat.anno, Periodo ~ series, fun.aggregate = sum) colnames(dat.anno) <- c("anno", "n.huelgas", "n.trabajadores", "n.jornadas") p <- ggplot(data = dat.anno, aes(x = anno, y = n.huelgas, group = rep(1, nrow(dat.anno)))) + geom_line() p <- p + geom_point(aes(size = n.jornadas)) p <- p + scale_x_discrete("año") + scale_y_continuous("número de huelgas") p ggsave("huelgas_por_anno.png") p <- ggplot(data = dat.anno, aes(x = anno, y = n.trabajadores/n.huelgas, group = rep(1, nrow(dat.anno)))) + geom_line() p <- p + scale_x_discrete("año") + scale_y_continuous("número de trabajadores por huelga") p ggsave("trabajadores_huelga_por_anno.png") p <- ggplot(data = dat.anno, aes(x = anno, y = n.jornadas /n.huelgas, group = rep(1, nrow(dat.anno)))) + geom_line() p <- p + scale_x_discrete("año") + scale_y_continuous("número de jornadas por huelga") p ggsave("jornadas_huelga_anno.png") para obtener, por un lado, el número de huelgas por mes desde enero de 1995 a noviembre de 2011: ...

29 de marzo de 2012 · Carlos J. Gil Bellosta

R y la distribución de Rayleigh

En la reunión de usuarios de R de Madrid de ayer, Carlos Ortega estudió la distribución en el tiempo del número de bugs que aparecen en el código de R en cada versión. Indicó que es plausible que sigan una distribución de Rayleigh, relativamente frecuente en ese tipo de contextos. E indicó que esta distribución, no tan conocida, tiene que ver (he olvidado lo que dijo exactamente) con dos normales independientes. Efectivamente, según la Wikipedia, la distribución de Rayleigh (de parámetro $\sigma$)admite la caracterización ...

23 de marzo de 2012 · Carlos J. Gil Bellosta