Sobre la peculiarísima implementación del modelo lineal en (pseudo-)scikit-learn

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.

2 comentarios sobre “Sobre la peculiarísima implementación del modelo lineal en (pseudo-)scikit-learn

  1. Kobayashi 17 julio, 2019 11:48

    No se si será por la versión de sklearn, pero no he conseguido replicar los resultados con la versión 0.21.2


    n = 1000
    X = np.random.rand(n, 2)
    X = np.column_stack((X, 1 * X[:,1]))
    Y = np.dot(X, np.array([1, 2, 1])) + 1 + np.random.randn(n) / 2
    reg = LinearRegression().fit(X, Y)
    reg.coef_

    #array([0.97399435, 1.47193515, 1.47193515])

    No obstante, sí es cierto lo de que internamente, el cálculo de los parámetros lo hace con un optimizador sobre la métrica del error en mínimos cuadrados, ya que diferentes ejecuciones devuelven diferentes resultados y no debería, ya que la regresión lineal es un problema de optimización que tiene solución exacta.

    Muy interesante el post.

    Saludos

  2. Carlos J. Gil Bellosta 17 julio, 2019 15:47

    Curioso que haya resultados muy distintos (cualitativamente) con versiones distintas. De todos modos, no sé concretamente a qué te refieres con que obtienes resultados distintos en ejecuciones distintas. Los resultados pueden bailar un poco si no fijas la semilla: la Y se genera de manera distinta (pero parecida) en cada ejecución del bloque.

    Por otro lado, con datos como los que he construido, la solución no es única: está en la región del espacio en que se cumple que a2 + a3 = 3 (donde ai es el coeficiente i-ésimo de la regresión).

Comenta

Your email address will not be published.

Puedes usar estas etiquetas y atributos de HTML: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong>

Este sitio usa Akismet para reducir el spam. Aprende cómo se procesan los datos de tus comentarios.