Para explorar como a LASSOregressão funciona, escrevi um pequeno pedaço de código que deve otimizar a LASSOregressão escolhendo o melhor parâmetro alfa.

Não consigo descobrir por que a LASSOregressão está me dando resultados tão instáveis ​​para o parâmetro alfa após a validação cruzada.

Aqui está o meu código 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()

Se você executar esse código repetidamente, ele fornecerá resultados totalmente diferentes para alfa:

>>> 
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 que o valor alfa não está convergindo corretamente? Eu sei que meus dados são sintéticos, mas a distribuição é a mesma. Além disso, a variação é muito pequena em x1e x2.

o que poderia estar causando isso tão instável?

A mesma coisa escrita em R fornece resultados diferentes - sempre retorna o valor mais alto possível para alfa como o "ideal_alpha".

Também escrevi isso em R, o que me dá uma resposta um pouco diferente, e não sei 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))

A saída do código R acima é:

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

De fato, não importa o que eu defina para a linha " alphas = logspace(-5, 6, 100)", sempre recebo o valor mais alto de alfa.

Eu acho que existem duas perguntas diferentes aqui:

1. Por que o valor alfa é tão instável para a versão escrita em Python?

2. Por que a versão escrita em R me dá um resultado diferente? (Sei que a logspacefunção é diferente Rpara python, mas a versão escrito em Rsempre me dá o maior valor de alphapara o valor alfa ideal, enquanto que a versão python não).

Seria ótimo saber essas coisas ...

Não conheço python muito bem, mas encontrei um problema com seu código R.

Você tem as 2 linhas:

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

Que soma os resíduos e depois os esquadrinha. Isso é muito diferente de agrupar os resíduos e depois somar (o que parece que o código python faz corretamente). Eu suspeitaria que isso possa ser uma grande parte da diferença entre o código R e o código python. Corrija o código R e execute-o novamente para ver se ele se comporta mais como o código python.

Eu também sugeriria que, em vez de salvar o "melhor" alfa e o mse correspondente, você armazene todos eles e plote o relacionamento. Pode ser que, para a sua configuração, exista uma região bastante plana, para que a diferença entre o mse em diferentes pontos não seja muito grande. Se for esse o caso, alterações muito pequenas nos dados (mesmo a ordem na validação cruzada) podem alterar qual ponto, dentre muitos essencialmente iguais, fornece o mínimo. Ter uma situação que resulte em uma região plana em torno do ideal geralmente leva ao que você está vendo, e o gráfico de todos os valores alfa com os valores mse correspondentes pode ser esclarecedor.

O sklearn tem um exemplo quase idêntico ao que você está tentando fazer aqui: http://scikit-learn.org/stable/auto_examples/exercises/plot_cv_diabetes.html

De fato, este exemplo mostra que você obtém resultados bastante variados para alfa para cada uma das três dobras realizadas nesse exemplo. Isso significa que você não pode confiar na seleção de alfa, porque é claramente altamente dependente da parte de seus dados que você está usando para treinar e selecionar alfa.

Não acho que você deva pensar na validação cruzada como algo que "convergirá" para lhe dar uma resposta perfeita. Na verdade, acho que conceitualmente é quase o oposto de convergir. Você está separando seus dados e, para cada dobra, está indo em uma 'direção separada'. O fato de você obter resultados diferentes, dependendo de como particionar seus dados de teste e treinamento, deve informar que a convergência em um resultado perfeito é impossível - e também não é desejável. A única maneira de obter um valor alfa consistente o tempo todo é se você usar todos os seus dados para treinamento. No entanto, se você fizer isso, obterá o melhor resultado de aprendizado, mas o pior resultado de validação.

A multicolinearidade x1e x2é o que torna o valor instável no código Python. A variação é tão pequena para as distribuições que geram essas variáveis, para que a variação dos coeficientes seja inflada. O fator de inflação de variação (VIF) pode ser calculado para ilustrar isso. Depois que a variação é aumentada de$\alpha$

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

....para....

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

então o valor estabiliza.$\alpha$

No Rentanto, o problema com o código ser diferente do código Python ainda é um mistério ...

Vou comentar sobre o código R:

Você está redefinindo variáveis ​​nos lugares errados, ou seja, as variáveis min_msedevem ser inicializadas como Inffora do forloop e optimal_alphadevem ser inicializadas como NULLlá. Isso se torna:

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)

A saída agora deve ser consistentemente os menores valores de alfa, porque há forte colinearidade nos preditores e a resposta é construída apenas a partir dos preditores disponíveis, ou seja, não há variável redundante que queremos que o LASSO coloque em zero, neste caso, não deseja realizar a regularização, ou seja, o menor alphadeve ser o melhor. Você pode ver o efeito do MSE aqui:

Observe que estou usando 50 alfas na mesma escala que você. Em torno do índice alfa 35, ambas as variáveis ​​são zeradas, o que significa que o modelo está sempre fazendo a mesma coisa e o mse fica estagnado.

Um problema melhor para estudar MSE, CV e LASSO

O problema acima não é muito interessante para o LASSO. O LASSO realiza a seleção de modelos, por isso queremos ver realmente escolher os parâmetros de interesse. É mais impressionante ver que o modelo está realmente escolhendo um alfa que realmente reduz o MSE, ou seja, nos dá melhores previsões ao jogar fora algumas variáveis. Aqui está um exemplo melhor, onde adiciono vários preditores redundantes.

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()

Então você apenas executa o loop for como no código acima! Note que eu coloquei o máximo do valor alphasde 1,5 para 6 a 1,5, apenas para ver o efeito no gráfico abaixo. Agora, o melhor alphavalor não é o mais baixo, mas você pode ver no gráfico que o MSE de validação cruzada está caindo e aumentando novamente no final. O ponto mais baixo desse gráfico corresponde ao índice alfa com o menor erro CV.