¿Por qué lazo no convergen en un parámetro de penalización?

Question

Para explorar cómo la LASSO de regresión obras, escribí un pequeño trozo de código que se debe optimizar LASSO de regresión escogiendo el mejor parámetro alfa.

No puedo entender por qué el LASSO de regresión me está dando tan inestable resultados para el parámetro alfa después de la validación cruzada.

Aquí está mi código de Python:

from sklearn.linear_model import Lasso
from sklearn.cross_validation import KFold
from matplotlib import pyplot as plt

# generate some sparse data to play with
import numpy as np
import pandas as pd 
from scipy.stats import norm
from scipy.stats import uniform

### generate your own data here

n = 1000

x1x2corr = 1.1
x1x3corr = 1.0
x1 = range(n) + norm.rvs(0, 1, n) + 50
x2 =  map(lambda aval: aval*x1x2corr, x1) + norm.rvs(0, 2, n) + 500
y = x1 + x2 #+ norm.rvs(0,10, n)

Xdf = pd.DataFrame()
Xdf['x1'] = x1
Xdf['x2'] = x2

X = Xdf.as_matrix()

# Split data in train set and test set
n_samples = X.shape[0]
X_train, y_train = X[:n_samples / 2], y[:n_samples / 2]
X_test, y_test = X[n_samples / 2:], y[n_samples / 2:]

kf = KFold(X_train.shape[0], n_folds = 10, )
alphas = np.logspace(-16, 8, num = 1000, base = 2)

e_alphas = list()
e_alphas_r = list()  # holds average r2 error
for alpha in alphas:
    lasso = Lasso(alpha=alpha, tol=0.004)
    err = list()
    err_2 = list()
    for tr_idx, tt_idx in kf:
        X_tr, X_tt = X_train[tr_idx], X_test[tt_idx]
        y_tr, y_tt = y_train[tr_idx], y_test[tt_idx]
        lasso.fit(X_tr, y_tr)
        y_hat = lasso.predict(X_tt)

        # returns the coefficient of determination (R^2 value)
        err_2.append(lasso.score(X_tt, y_tt))

        # returns MSE
        err.append(np.average((y_hat - y_tt)**2))
    e_alphas.append(np.average(err))
    e_alphas_r.append(np.average(err_2))

## print out the alpha that gives the minimum error
print 'the minimum value of error is ', e_alphas[e_alphas.index(min(e_alphas))]
print ' the minimizer is ',  alphas[e_alphas.index(min(e_alphas))]

##  <<< plotting alphas against error >>>

plt.figsize = (15, 15)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(alphas, e_alphas, 'b-')
ax.plot(alphas, e_alphas_r, 'g--')
ax.set_ylim(min(e_alphas),max(e_alphas))
ax.set_xlim(min(alphas),max(alphas))
ax.set_xlabel("alpha")
plt.show()

Si ejecuta este código en repetidas ocasiones, se da totalmente diferentes resultados para alpha:

>>> 
the minimum value of error is  3.99254192539
 the minimizer is  1.52587890625e-05
>>> ================================ RESTART ================================
>>> 
the minimum value of error is  4.07412455842
 the minimizer is  6.45622425334
>>> ================================ RESTART ================================
>>> 
the minimum value of error is  4.25898253597
 the minimizer is  1.52587890625e-05
>>> ================================ RESTART ================================
>>> 
the minimum value of error is  3.79392968781
 the minimizer is  28.8971008254
>>> 

¿Por qué es el valor de alfa no converge correctamente? Sé que mis datos es sintético, pero la distribución es la misma. También, la variación es muy pequeña en x1 y x2.

lo que podría estar causando este tan inestables?

Lo mismo escrito en R da resultados diferentes - siempre devuelve el valor más alto posible para el alfa como el "optimal_alpha".

También escribí esto en R, lo que me da un poco diferente de respuesta, que no sé por qué?

library(glmnet)
library(lars)
library(pracma)

set.seed(1)
k = 2 # number of features selected 

n = 1000

x1x2corr = 1.1
x1 = seq(n) + rnorm(n, 0, 1) + 50
x2 =  x1*x1x2corr + rnorm(n, 0, 2) + 500
y = x1 + x2 

filter_out_label <- function(col) {col!="y"}

alphas = logspace(-5, 6, 100)

for (alpha in alphas){
  k = 10
  optimal_alpha = NULL
  folds <- cut(seq(1, nrow(df)), breaks=k, labels=FALSE)
  total_mse = 0
  min_mse = 10000000
  for(i in 1:k){
    # Segement your data by fold using the which() function
    testIndexes <- which(folds==i, arr.ind=TRUE)
    testData <- df[testIndexes, ]
    trainData <- df[-testIndexes, ]

    fit <- lars(as.matrix(trainData[Filter(filter_out_label, names(df))]),
                trainData$y,
                type="lasso")
    # predict
    y_preds <- predict(fit, as.matrix(testData[Filter(filter_out_label, names(df))]),
                       s=alpha, type="fit", mode="lambda")$fit # default mode="step"

    y_true = testData$y
    residuals = (y_true - y_preds)
    mse=sum(residuals^2)
    total_mse = total_mse + mse
  }
  if (total_mse < min_mse){
    min_mse = total_mse
    optimal_alpha = alpha
  }
}

print(paste("the optimal alpha is ", optimal_alpha))

La salida de la R código anterior es:

> source('~.....')
[1] "the optimal alpha is  1e+06"

De hecho, no importa lo que me puse para la línea "alphas = logspace(-5, 6, 100)", siempre vuelvo el valor más alto de alfa.

Supongo que en realidad hay 2 diferentes preguntas aquí :

1. ¿Por qué es el valor de alfa de forma inestable, para la versión escrita en Python?

2. ¿Por qué la versión escrita en R me da un resultado diferente? (Me doy cuenta de que el logspace función es diferente de R a python, pero la versión escrita en R siempre me da el mayor valor de alpha , para el óptimo valor de alfa, mientras que la versión de python no).

Sería bueno saber estas cosas...

Answer 1

No sé python muy bien, pero me encontré un problema con tu código R.

Usted tiene 2 líneas:

residuals = sum(y_true - y_preds)
mse=residuals^2

Que las sumas de los residuos, a continuación, plazas de ellos. Esto es muy diferente de el cuadrado de los residuos, sumando después de ellos (el que aparece que el código python hace correctamente). Yo sospecho que esto puede ser una gran parte de la diferencia entre el R código y el código de python. Revisión de la R código y ejecutarlo de nuevo para ver si se comporta más como el código de python.

También me gustaría sugerir que en lugar de salvar el "mejor" el alfa y la correspondiente mse que almacenar todos ellos y de la trama de la relación. Podría ser que para su instalación no es una región que es bastante plana, con lo que la diferencia entre el mse en diferentes puntos no es muy grande. Si este es el caso, entonces cambios muy menores a los datos (incluso el orden en la validación cruzada) puede cambiar el punto de que, entre los muchos que son esencialmente el mismo, da el mínimo. Tener una situación que se traduce en una región plana alrededor de la óptima frecuencia conduce a lo que están viendo y la trama de todos los valores alfa con los correspondientes valores de mse podría ser esclarecedor.

Answer 2

sklearn tiene un ejemplo que es casi idéntico a lo que estamos tratando de hacer aquí: http://scikit-learn.org/stable/auto_examples/exercises/plot_cv_diabetes.html

De hecho, este ejemplo muestra que usted consiga salvajemente diferentes resultados para alfa para cada uno de los tres pliegues hecho en el ejemplo. Esto significa que usted no puede confiar en la selección de alfa porque claramente es altamente dependiente de lo que una porción de los datos que está utilizando para entrenar y seleccione alfa.

No creo que usted debe pensar en la validación cruzada como algo que va a "convergencia", para darle una respuesta perfecta. En realidad, creo que conceptualmente es casi el opuesto de la convergencia. Se están separando de sus datos y para cada una de las veces que se va en un 'separado de dirección". El hecho de que usted obtenga resultados diferentes dependiendo de cómo particionar su entrenamiento y de prueba de datos debe de decirle que convergen en un resultado perfecto es imposible, y además no es deseable. La única manera en la que podría obtener una constante de valor alfa todo el tiempo es que si usted fuera a utilizar todos sus datos para el entrenamiento. Sin embargo, si usted fuera a hacer esto, usted tendrá la mejor forma de aprender resultado, pero el peor resultado de la validación.

Answer 3

El multi-colinealidad en x1 y x2 es lo que hace que el $\alpha$ valor inestable en el código de Python. La variación es tan pequeña para las distribuciones que generan estas variables, de modo que la varianza de los coeficientes se infla. El factor de inflación de la varianza (VIF) puede ser calculada para ilustrar esto. Después de que la varianza es mayor de

x1 = range(n) + norm.rvs(0, 1, n) + 50
x2 =  map(lambda aval: aval*x1x2corr, x1) + norm.rvs(0, 2, n) + 500

....a....

x1 = range(n) + norm.rvs(0, 100, n) + 50
x2 =  map(lambda aval: aval*x1x2corr, x1) + norm.rvs(0, 200, n) + 500

a continuación, el $\alpha$ valor se estabiliza.

El problema con la R código diferente del código de Python es todavía un misterio, sin embargo...

Answer 4

Voy a comentar sobre el código R:

Restablece las variables en los lugares equivocados, es decir, las variables min_mse debe ser inicializado como Inf fuera de la for de bucle y optimal_alpha debe ser inicializado como NULL no. Esto se convierte en:

library(glmnet)
library(lars)
library(pracma)

set.seed(1)
k = 2 # number of features selected 

n = 100

x1x2corr = 1.1
x1 = seq(n) + rnorm(n, 0, 1) + 50
x2 =  x1*x1x2corr + rnorm(n, 0, 2) + 500
y = x1 + x2 +rnorm(n,0,0.5)
df = data.frame(x1 = x1, x2 = x2, y = y)
filter_out_label <- function(col) {col!="y"}

alphas = logspace(-5, 6, 50)

###
# INITIALIZE here before loop
###
min_mse = Inf
optimal_alpha = NULL
# Let's store the mse values for good measure
my_mse = c()

for (alpha in alphas){
  k = 10
  folds <- cut(seq(1, nrow(df)), breaks=k, labels=FALSE)
  # DO NOT INITIALIZE min_mse and optimal_alpha here, 
  # then you cannot find them...
  total_mse = 0
  for(i in 1:k){
    # Segement your data by fold using the which() function
    testIndexes <- which(folds==i, arr.ind=TRUE)
    testData <- df[testIndexes, ]
    trainData <- df[-testIndexes, ]

    fit <- lars(as.matrix(trainData[Filter(filter_out_label, names(df))]),
                trainData$y,
                type="lasso")
    # predict
    y_preds <- predict(fit, as.matrix(testData[Filter(filter_out_label,
                       names(df))]),
                       s=alpha, type="fit", mode="lambda")$fit 

    y_true = testData$y
    residuals = (y_true - y_preds)
    mse=sum(residuals^2)
    total_mse = total_mse + mse
  }
  # Let's store the MSE to see the effect
  my_mse <- c(my_mse, total_mse)
  if (total_mse < min_mse){
    min_mse = total_mse
    optimal_alpha = alpha
    # Let's observe the output
    print(min_mse)
  }
}

print(paste("the optimal alpha is ", optimal_alpha))
# Plot the effect of MSE with varying alphas
plot(my_mse)

La salida debe ahora ser siempre el más pequeño de los valores de alfa, porque hay una fuerte colinealidad en los predictores y la respuesta sólo es construido a partir de la disposición de los predictores, es decir, no hay redundante de la variable que queremos LAZO para poner a cero, en este caso no queremos realizar la regularización, es decir, el más pequeño alpha debe ser el mejor. Usted puede ver el efecto de MSE aquí:

Tenga en cuenta que estoy usando 50 alfas en la misma escala. Alrededor de alfa indexado 35 ambas variables se cerró de golpe a cero, lo que significa que el modelo es hacer siempre lo mismo y el mse se estanca.

Un problema mejor para el estudio de MSE, CV y el LAZO

El problema anterior no es muy interesante para el LAZO. El LAZO realiza la selección del modelo, por lo que queremos ver es realmente elegir los parámetros de interés. Es más impresionante ver que el modelo es en realidad la selección de un alfa que, en realidad, reduce el MSE, es decir, que nos da mejores predicciones de tirar de algunas variables. Aquí hay un ejemplo mejor, donde puedo añadir un montón de redundante predictores.

set.seed(1)
k = 100 # number of features selected 

n = 100

x1x2corr = 1.1
x1 = seq(n) + rnorm(n, 0, 1) + 50
x2 =  x1*x1x2corr + rnorm(n, 0, 2) + 500
# Rest of the variables are just noise
x3 = matrix(rnorm(k-2,0,(k-2)*n),n,k-2)
y = x1 + x2 +rnorm(n,0,0.5)
df = data.frame(x1 = x1, x2 = x2, y = y)
df <- cbind(df,x3)
filter_out_label <- function(col) {col!="y"}

alphas = logspace(-5, 1.5, 100)
min_mse = Inf
optimal_alpha = NULL
my_mse = c()

Entonces usted acaba de ejecutar el bucle for como en el código de arriba! Tenga en cuenta que pongo el máximo de la alphas 1,5 de 6, solo para ver el efecto en el gráfico siguiente. Ahora el mejor alpha valor no es el más bajo, pero se puede ver en el gráfico que el de validación cruzada MSE es tomar una gota y picos de nuevo en la final. El punto más bajo en ese gráfico, se corresponde con el índice alfa con el menor CV-error.