R

Me escriben pidiendo consejo sobre cómo leer datos contenidos en (una serie larga de) ficheros en formatos .dbf, .xlsx (con un formato extraño) y .pdf.

.dbf

No tengo ni curiosidad por averiguar de dónde proceden. Simplemente,

library(foreign)
res <-read.dbf("R0010.DBF")

funciona de maravilla.

.xlsx

Estos sí que sé de dónde vienen (y me guardo la opinión). El problema aquí no era leer directamente tablas contenidas en hojas sino ir extrayendo celdas y rangos de hojas. Así que:

Se ve que hay arqueólogos bayesianos. Un problema con el que se encuentran es que tropiezan con cacharros antiguos y quieren estimar su antigüedad.

Así que prueban distintos métodos (¿químicos?), cada uno de los cuales con su precisión, y acaban recopilando una serie de estimaciones y errores. Obviamente, tienen que combinarlas de alguna manera.

El modelo más simple es

$$ M_i \sim N(\mu, \sigma_i)$$

donde $latex \mu$ es la antigüedad (desconocida) del artefacto y los $latex \sigma_i$ son las varianzas distintas de los distintos métodos de medida, que arrojan las estimaciones $latex M_i$.

Envalentonado por el comentario de Iñaki Úcar a mi entrada del otro día, que me remitía a este artículo, decidí rizar el rizo y crear intervalos de confianza no ya discontinuos sino con otra propiedad topológica imposible: homeomorfos con un toro.

Y aquí está:

El modelo, el código y demás,

library(rstan)
library(ggplot2)

n <- 100

a1 <- 1
a2 <- 1
sigma <- 0.4

datos <- data.frame(x1 = rnorm(n, 2, 0.1),
                    x2 = rnorm(n, 2, 0.1))

datos$y <- a1^datos$x1 + a2^datos$x2 + rnorm(n, 0, sigma)

codigo <- "
data {
  int<lower=1> N;
  real y[N];
  real x1[N];
  real x2[N];
}

parameters {
  real<lower=-3, upper="3"> a1;
  real<lower=-3, upper="3"> a2;
  real<lower=0, upper="3"> sigma;
}

model {
  for (n in 1:N)
    y[n] ~ normal(fabs(a1)^x1[n] +
      fabs(a2)^x2[n], sigma);
}"

fit <- stan(model_code = codigo,
    data = list(N = length(datos$y), y = datos$y,
                x1 = datos$x1, x2 = datos$x2),
    iter=40000, warmup=2000,
    chains=1, thin=10)

res <- as.data.frame(fit)

ggplot(res, aes(x = a1, y = a2)) + geom_point(alpha = 0.1)

De nuevo, no son intervalos propiamente dichos, lo convengo. Pero son configuraciones más fieles al espíritu de lo que un intervalo de confianza es y representa que su(s) letra(s) I N T E R V A L O.

    curve(-sqrt(x^2 + 1), -5, 5)

pinta una rama de hipérbola,

que, una vez exponenciada, i.e.,

    curve(exp(-sqrt(x^2 + 1)), -5, 5)

da

Es decir, una curva algo menos esbelta que la normal pero que bien podemos dividir por su integral para obtener la llamada distribución hiperbólica.

Tres notas sobre ella:

• Tiene una historia curiosa. Fue considerada por Ralph Bagnold al estudiar la forma de las dunas y la sedimentación de la arena arrastrada por el viento. El logaritmo de sus curvas, se ve, tenía forma de hipérbola.
• Lo cual os proporciona un exótico contraejemplo al argumento habitual sobre la naturaleza omniatractora de la normal.
• La distribución hiperbólica (y sus extensiones) están disponibles en el paquete ghyp, motivado por aplicaciones financieras, como siempre. Esa gente es adicta a distribuciones con colas gruesas. Aunque para lo que les valen luego…

Continúo con esto que concluí con una discusión que me negué a resolver sobre la geometría de los errores.

Que es la manera de entender que los problemas directos e inversos no son exactamente el mismo. Digamos que no es una medida invariante frente a reflexiones del plano (que es lo que hacemos realmente al considerar el modelo inverso).

¿Pero y si medimos la distancia (ortogonal) entre los puntos $latex (x,y)$ y la curva $latex y = f(x)$ (o, equivalentemente, $latex x = f^{-1}(x)$)? Entonces daría (o debería dar) lo mismo.

Inspirado por esto he tratado de contrastar una hipótesis en otro contexto.

Las cosas, o se hacen bien, o no se hacen. Como mi análisis se ha complicado con casos y casitos particulares, aunque siga pensándo cierta (en caso de tener que apostar, como priori, claro) la hipótesis de partida, abandono su búsqueda.

Como subproducto, esto:

library(xml2)
library(stringr)
library(plyr)
library(lubridate)

periodos <- expand.grid(anno = 2010:2017, mes = 1:12)
periodos$ind <- periodos$anno * 100 + periodos$mes
periodos <- periodos[periodos$ind < 201711,]
periodos <- paste(periodos$anno,
  str_pad(periodos$mes, 2, pad = "0"), sep = "_")

raw <- lapply(periodos, function(x){
  url <- paste0("http://www.eldiario.es/sitemap_contents_", x, ".xml")
  print(url)
  as_list(read_xml(url))
})

#df <- lapply(raw, function(y)
  ldply(y, function(x) as.data.frame(t(unlist(x)))))

res <- lapply(raw, unlist)
res <- lapply(res, function(x) t(matrix(x, 3, length(x) / 3)))
res <- data.frame(url = res[,1],
  time = res[,2], stringsAsFactors = FALSE)

res$time <- gsub("\\+.*", "", res$time)
res$time <- strptime(res$time,
  "%Y-%m-%dT%H:%M:%S")

res$titular <- gsub("_0_[0-9]*.html", "", res$url)
res$titular <- gsub(".*/", "", res$titular)
res$titular <- tolower(res$titular)

res$year <- year(res$time)
res$month <- month(res$time)

Igual le sirve a alguien para analizar palabras clave en titulares de ese u otro medio, su evolución por mes, etc.

Las malandanzas de Circiter la han conducido al siguiente entuerto: estimar $latex \alpha$ donde

$$ y = f_\alpha(x) + \epsilon$$

y $latex f_\alpha$ es una función no lineal horrible. Sin embargo, $latex f^{-1}_\alpha$ es mucho más manejable y podría plantearse el modelo

$$ x = f^{-1}_\alpha(y) + \epsilon$$

(donde este nuevo $latex \epsilon$ no coincide con el anterior: piénsese en el método delta y léase la nota final).

Un ejemplo. Que arranca con unos datos autoexplicativos:

La gente de rOpenSci hace cosas a las que merece la pena atento. Tanto por los objetivos como por medios y las formas. Recomiendo seguir sus últimas publicaciones acerca de la profesionalización del proceso de desarrollo de código.

Llevo unos meses jugando con una idea inspirada por rOpenSci: crear un respositorio y un consorcio más o menos formal que desarrolle, mantenga y mejore herramientas (en R) de interés para el procesamiento y análisis de datos ya no científicos sino españoles. Hablo, obviamente, de INE (EPA, EPF, censo, padrón,…), CIS (barómetros, etc.), IGN (Siane,…), encuestas electorales, etc.

Una función no debería cambiar nada de cuanto la rodea. Debería devolver algo y ya. Se acepta barco como animal acuático cuando hay funciones que escriben en logs, guardan datos en disco o crean gráficos.

R deja que los usuarios se disparen en el pie permitiendo hacer cosas tan peligrosas como:

a <- new.env()

a$1     # error

foo <- function(){
  a$a <- 1
}

foo()
a$a
# [1] 1

De la misma manera, si le enseñas un cuchillo a una vieja, es posible que te dé su bolso con todo lo que contiene. Pero eso no significa que debas usar los cuchillos para tales fines.

El título, no el de esta entrada sino el de A Comparison of Programming Languages in Economics, es una sinécdoque confusa.

Que nadie busque en él consejo sobre qué lenguaje estudiar si le interesa el mundo de la economía (en general). O fuera de ella (también en general).

Encontrará más bien la implementación de la solución a un único problema dentro de los muchos que supongo comprende esa disciplina. Uno, además, con el que no he visto (en persona) a economista alguno ganarse el pan ni en la academia ni fuera de ella.