Python

Hoy me he retrasado en escribir por haber estado probando (y estresando, como hay quien dice), software para resolver problemas de programación lineal. En total, nada, unos diez millones de variables unos treinta millones de restricciones.

Nota: es un problema LP puro, nada de enteros, nada de pérdidas no lineales, etc.

• Primera opción: Python + PuLP + CBC (de COIN-OR), que es el optimizador por defecto de PuLP. Rendimiento aceptable para el tipo de uso que se le acabaría dando. Se ha convertido en el benchmark.
• Segunda opción: Python + OR-Tools (de Google), y en particular, Glop. Un tanto decepcionante: aunque ne términos de velocidad no es apreciablemente inferior a CBC, en muchos casos desistía y no encontraba ninguna solución.

Este tipo de problemas y yo nos reencontramos indefectiblemente cada cinco años. Así que, de una vez a otra, se me ha olvidado casi todo. De modo que si alguien tiene el asunto más fresco y le da rabia que algún diletante como opte por soluciones subóptimas y/o viejunas y esté entre asombrado e indignado de que ignore el último grito de la cosa, tiene la posibilidad de enmendarme a mí y enseñarnos, de paso, a todos, en los comentarios.

Así vendría a traducirse el título de este artículo, que trata de taxonomizar y sistematizar una serie de técnicas muy recientes para explicar modelos de caja negra.

Tal vez no acabe siendo la manera pero, sin duda, acabará habiendo una.

Así arranca este artículo, que presenta una extensión de XGBoost para predicciones probabilísticas. Es decir, un paquete que promete no solo una estimación del valor central de la predicción sino de su distribución.

La versión equivalente de lo anterior en el mundo de los random forests está descrito aquí, disponible aquí y mucho me temo que muy pronto voy a poder contar por aquí si está a la altura de las expectativas.

Se acaba de publicar en GitHub el/nuestro Curso de python básico orientado al análisis de datos.

Digo nuestro un tanto impropiamente: casi todo el material es de Luz Frías, mi socia en Circiter. Mía hay alguna cosa suelta.

Como como minicoautor soy el comentarista menos creíble del contenido, lo dejo al juicio de cada cual. Y, por supuesto, se agradecen correcciones, comentarios, cañas y fusilamientos (con la debida caballerosidad, por supuesto, en lo de las atribuciones).

Este martes 17 de diciembre hablaré durante una hora sobre (cierto tipo de) big data y modelos adecuados para modelizarlos en el CodingClub de la Universidad Carlos III.

• El contenido de la charla, entiendo, se publicará también después en el blog del CodingClub.
• Los detalles (sitio, hora, etc.) están en el enlace indicado más arriba.
• Obviamente, agradezco a los organizadores del CodingClub por haberme invitado. Espero no estar arrepentido el martes por la tarde de lo siguiente: es el ciclo de charlas sobre cosas relacionadas con datos más seria y mejor organizada que conozco.

Y con eso, prácticamente, cierro el 2019 para casi todos los efectos. En 2020, más.

Pensé que ya había escrito sobre el asunto porque tropecé con él en un proyecto hace un tiempo. Pero mi menoria se había confundido con otra entrada, Sobre la peculiarisima implementacion del modelo lineal en (pseudo-)Scikit-learn, donde se discute, precisamente, un problema similar si se lo mira de cierta manera o diametralmente opuesto si se ve con otra perspectiva.

Allí el problema era que Scikit-learn gestionaba muy sui generis el insidioso problema de la colinealidad. Precisamente, porque utiliza un optimizador ad hoc y no estándar para ajustar el modelo lineal.

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.