Python

Leyendo sobre si dizque PyTorch le siega la hierba debajo de los pies a TensorFlow, averigué la existencia de Pyro.

Pyro se autopresenta como Deep Universal Probabilistic Programming, pero aplicando métodos porfirianos (ya sabéis: género próximo y diferencia específica), es, o pretende ser, Stan en Python y a escala.

Aquí van mis dos primeras impresiones, basadas en una inspección superficial de los tutoriales.

En primer lugar, aunque Pyro permite usar (distintas versiones de) MCMC, parece que su especialidad es la inferencia variacional estocástica. Que parece funcionar de la siguiente manera. En el MCMC tradicional uno obtiene una muestra de la distribución (a posteriori, para los amigos) de los parámetros de interés. Eso es todo: vectores de puntos. En la inferencia variacional estocástica, uno preespecifica la forma paramétrica de la posteriori y el algoritmo calcula sus parámetros a partir de los valores simulados. Por ejemplo, uno va y dice: me da que la distribución del término independiente de mi regresión lineal va a ser normal. Entonces, Pyro responde: si es normal, la mejor media y desviación estándar que encuentro son tal y cual.

Es una pregunta que surge reiteradamente. Por ejemplo, cuando se compara un clúster con el resto de la población y uno busca las variables que mejor lo caracterizan. Y crear gráficos como

(extraído de aquí) donde las variables están ordenadas de acuerdo con su poder discriminador.

Mi técnica favorita para crear tales indicadores es la EMD (earth mover’s distance) y/o sus generalizaciones, muy bien descritas en Optimal Transport and Wasserstein Distance y disponibles en R y Python.

Los modelos mixtos en Python son un bien público.

El sector privado no produce suficientes bienes públicos (con excepciones tan notables como las búsquedas en Google o las páginas aún sin paywall de los periódicos). El sector público y los impuestos que lo financian argumenta la conveniencia de su propia existencia en términos de esa provisión de bienes públicos que dizque realiza.

Pero ese subsector del sector público que debería implementar los modelos mixtos en Python se dedica a otra cosa. A cosas como estas.

Si ejecutas

import numpy as np
from sklearn.linear_model import LinearRegression

n = 1000
X = np.random.rand(n, 2)

Y = np.dot(X, np.array([1, 2])) + 1 + np.random.randn(n) / 2
reg = LinearRegression().fit(X, Y)

reg.intercept_
reg.coef_

se obtiene más o menos lo esperado. Pero si añades una columna linealmente dependiente,

X = np.column_stack((X, 1 * X[:,1]))

ocurren cosas de la más calamitosa especie:

Y = np.dot(X, np.array([1, 2, 1])) + 1 + np.random.randn(n) / 2
reg = LinearRegression().fit(X, Y)
reg.coef_
# array([ 9.89633991e-01, -1.63740303e+14,  1.63740303e+14])

Comentarios:

• Diríase que la implementación del modelo lineal en scikit-learn no es la que se estudia por doquier (la prima, la inversa, etc.); sospecho que convierte el error cuadrático en una función que depende de los coeficientes y se la pasan a un optimizador (más o menos) genérico.
• Supongo que la implementación actual pasa todos las pruebas unitarias.
• Y sospecho, además, que las pruebas unitarias no las ha planteado un estadístico.
• Así que tal vez los de scikit-learn no saben que tienen problemas de colinealidad; y si alguien se lo ha comentado, igual no han comprendido el issue.
• Para la predicción, el problema arriba apuntado no es tal. Aun con coeficientes desaforados y siempre que no haya problemas de precisión numérica, tanto da hacer las cosas como todo el mundo o implementando ocurrencias como la anterior.
• Pero para todo lo demás (p.e., impacto de variables, etc.), la implementación es de traca y no vale para maldita de Dios la cosa.
• Aunque a estas alturas de siglo, ¿quién en su sano juicio usa el modelo lineal básico?
• Además, en la práctica, no hay problemas de colinealidad (ni aproximada o, mucho menos, exacta). Hummm…
• ¿O sí? Mi colega Cañadas ha escrito una entrada en su blog sobre la codificación de variables categóricas donde denuncia cómo en Python las funciones más habituales crean por defecto columnas linealmente dependientes por defecto (por no omitir el primer nivel). O sea, en Python, si no andas con cuidado, acabas con la suela llena de kk de perro.

Resumen:

• He decidido usar RStudio como IDE para Python. RStudio no es el mejor IDE para desarrollar, pero es incomparablemente mejor que cualquier otro IDE para explorar, etc. Funciona muy bien y solo puede mejorar.
• He decidido pasar de Jupyter. Los notebooks valen para lo que valen, pero no para lo que hago. En caso de necesidad, uso Rmarkdown con bloques de Python. De nuevo, funcionan muy bien y solo pueden mejorar.
• Finalmente, he decidido pasar de Anaconda. Tiene incompatibilidades con RStudio. Particularmente, cuando los módulos de Python tratan de cargar shared libraries. Los módulos de Anaconda tienen el vicio de buscarlos dentro del directorio de instalación, pero al lanzar el intérprete de Python a través de reticulate, en Linux parece que los busca en el sistema (por debajo de /usr/lib y similares). Y todo se rompe mucho. Mucho y muy, muy feo.

Así que uso los Python (3.7 cuando puedo, otras versiones cuando me obligan) del sistema. Pero la instalación del sistema es mínima. He creado varios environments ad hoc (y dentro de un directorio ad hoc para ellos) y obigo a reticulate a usarlos (vía use_virtualenv()) según conveniencia. En ellos tengo todas las dependencias (de numpy para arriba).

Era casi todavía el siglo XX cuando yo, desesperado por hacer cosas que consideraba normales y que SAS no me permitía, pregunté a un profesor por algo como C pero para estadística. Y el profesor me contó que conocía a alguien que conocía a alguien que conocía a alguien que usaba una cosa nueva que se llamaba R y que podía servirme.

Fue amor a primera vista, pero esa es otra historia. La relevante aquí es que volví a hablar con aquel profesor para agradecerle el consejo y, de paso, le pregunté que por qué no lo usaba él. Me contestó que porque en R no había modelos mixtos (aunque nlme es anterior, del 99; ¡a saber en qué estado se encontraba entonces!).

Participé el otro día en una cena con gente friqui. Constaté con cierto desasosiego cómo han virado los sujetos pasivos de nuestra indignación profesional a lo largo de los años.

Antaño, fueron los viejos que seguían apegados a la paleoinformática. Hogaño, los primíparos que usan Python y desdeñan R.

Tengo sentimientos encontrados y no sé qué más añadir.

DataSynthesizer (véase también el correspondiente artículo) es un programa en Python que:

1. Toma una tabla de datos (microdatos, de hecho) que contiene información confidencial.
2. Genera otra aleatoria pero que conserva (¿los conservará?) la estructura básica de la información subyacente (conteos, correlaciones, etc.).

Está pensado para poder realizar el análisis estadístico de (determinados) datos sin verlos propiamente.

Particularmente interesante es el algoritmo para preservar la correlación entre columnas.

[Nota: he aprovechado la entrada para acuñar el neologismo microdatado para referirme a quien figura en un fichero de microdatos.]

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.

Hoy he visto

aquí y he escrito

m <- matrix(c(74, 15, 10, 1, 11, 50, 38, 1,
            5, 4, 90, 1, 17, 4, 19, 60),
            4, 4, byrow = TRUE)
m <- m / 100

luego

m %*% m %*% m %*% m %*% m %*% m %*% m %*% m %*% m %*% m %*% m %*% m %*% m%*% m%*% m%*% m%*% m%*% m%*% m %*% m %*% m %*% m %*% m %*% m %*% m %*% m %*% m %*% m %*% m %*% m %*% m%*% m%*% m%*% m%*% m%*% m%*% m
#          [,1]      [,2]      [,3]       [,4]
#[1,] 0.1926676 0.1133218 0.6696203 0.02439024
#[2,] 0.1926647 0.1133206 0.6696245 0.02439024
#[3,] 0.1926638 0.1133202 0.6696258 0.02439024
#[4,] 0.1926675 0.1133218 0.6696205 0.02439025

y finalmente