Um herauszufinden, wie die LASSO
Regression funktioniert, habe ich einen kleinen Code geschrieben, der die LASSO
Regression durch Auswahl des besten Alpha-Parameters optimieren soll .
Ich kann nicht herausfinden, warum die LASSO
Regression nach der Kreuzvalidierung so instabile Ergebnisse für den Alpha-Parameter liefert.
Hier ist mein Python-Code:
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()
Wenn Sie diesen Code wiederholt ausführen, ergeben sich für Alpha völlig unterschiedliche Ergebnisse:
>>>
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
>>>
Warum konvergiert der Alpha-Wert nicht richtig? Ich weiß, dass meine Daten synthetisch sind, aber die Verteilung ist dieselbe. Auch die Variation ist in x1
und sehr gering x2
.
Was könnte dazu führen, dass dies so instabil ist?
Dasselbe, was in R geschrieben ist, führt zu unterschiedlichen Ergebnissen - es gibt immer den höchstmöglichen Wert für Alpha als "optimal_alpha" zurück.
Ich habe das auch in R geschrieben, was mir eine etwas andere Antwort gibt, die ich nicht weiß warum?
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))
Die Ausgabe des obigen R-Codes lautet:
> source('~.....')
[1] "the optimal alpha is 1e+06"
Unabhängig davon, was ich für die Zeile " alphas = logspace(-5, 6, 100)
" eingestellt habe, erhalte ich immer den höchsten Wert für Alpha zurück.
Ich denke, hier gibt es tatsächlich zwei verschiedene Fragen:
Warum ist der Alpha-Wert für die in Python geschriebene Version so instabil?
Warum gibt mir die in R geschriebene Version ein anderes Ergebnis? (Mir ist klar, dass sich die
logspace
Funktion vonR
bis unterscheidetpython
, aber die eingeschriebene VersionR
gibt mir immer den größten Wertalpha
für den optimalen Alpha-Wert, während die Python-Version dies nicht tut).
Es wäre toll, diese Dinge zu wissen ...
fit_intercept
Parameter beim Erstellen des Lasso-Modells zu verwenden.Antworten:
Ich kenne Python nicht sehr gut, aber ich habe ein Problem mit Ihrem R-Code gefunden.
Sie haben die 2 Zeilen:
Was die Residuen summiert, quadriert sie dann. Dies unterscheidet sich stark vom Quadrieren der Residuen und anschließenden Summieren (was anscheinend der Python-Code korrekt macht). Ich würde vermuten, dass dies ein großer Teil des Unterschieds zwischen dem R-Code und dem Python-Code ist. Korrigieren Sie den R-Code und führen Sie ihn erneut aus, um festzustellen, ob er sich eher wie der Python-Code verhält.
Ich würde auch vorschlagen, dass Sie nicht nur das "beste" Alpha und das entsprechende mse speichern, sondern alle speichern und die Beziehung zeichnen. Es kann sein, dass es für Ihr Setup eine Region gibt, die ziemlich flach ist, so dass der Unterschied zwischen der mse an verschiedenen Punkten nicht sehr groß ist. Wenn dies der Fall ist, können sehr geringfügige Änderungen an den Daten (sogar die Reihenfolge in der Kreuzvalidierung) ändern, welcher Punkt unter vielen, die im Wesentlichen gleich sind, das Minimum ergibt. Eine Situation, die zu einem flachen Bereich um das Optimum führt, führt häufig zu dem, was Sie sehen, und die Darstellung aller Alpha-Werte mit den entsprechenden mse-Werten kann aufschlussreich sein.
quelle
sklearn
Version über eine integrierte Kreuzvalidierung, wie @JennyLu hervorhob, sodass ein etwas anderer Fehler auftritt.sklearn hat ein Beispiel, das fast identisch mit dem ist, was Sie hier versuchen: http://scikit-learn.org/stable/auto_examples/exercises/plot_cv_diabetes.html
In der Tat zeigt dieses Beispiel, dass Sie für jede der drei in diesem Beispiel durchgeführten Falten sehr unterschiedliche Ergebnisse für Alpha erhalten. Dies bedeutet, dass Sie der Auswahl von Alpha nicht vertrauen können, da dies eindeutig stark davon abhängt, welchen Teil Ihrer Daten Sie zum Trainieren und Auswählen von Alpha verwenden.
Ich denke nicht, dass Sie Kreuzvalidierung als etwas betrachten sollten, das "konvergiert", um Ihnen eine perfekte Antwort zu geben. Eigentlich denke ich, dass es konzeptionell fast das Gegenteil von Konvergenz ist. Sie trennen Ihre Daten und gehen für jede Falte in eine "separate Richtung". Die Tatsache, dass Sie je nach Partitionierung Ihrer Test- und Trainingsdaten unterschiedliche Ergebnisse erhalten, sollte Ihnen sagen, dass die Konvergenz zu einem perfekten Ergebnis unmöglich - und auch nicht wünschenswert ist. Der einzige Weg, um jederzeit einen konsistenten Alpha-Wert zu erhalten, besteht darin, alle Ihre Daten für das Training zu verwenden. Wenn Sie dies jedoch tun würden, würden Sie das beste Lernergebnis, aber das schlechteste Validierungsergebnis erhalten.
quelle
sklearn
Beispiel, das Sie auf plot_cv_diabetes zitiert haben, hat es so wenige Datenpunkte (150), dass ich nicht davon überzeugt wäre, dass diek folds cross-validation
und (so wie ich es im obigen Beispiel gemacht habe) der Wert vonk folds cross-validation
im Wesentlichen, den bedingten Mittelwert des von Ihnen geschätzten Parameters (den Schätzer) zu berechnen. Also denke ich nicht den Wert vonDie Multikollinearität inα Wert instabil im Python-Code. Die Varianz ist für die Verteilungen, die diese Variablen erzeugen, so gering, dass die Varianz der Koeffizienten aufgeblasen wird. Der Varianzinflationsfaktor (VIF) könnte berechnet werden, um dies zu veranschaulichen. Nachdem die Varianz von erhöht ist
x1
undx2
macht das aus....zu....
dann ist dieα Wert stabilisiert sich.
Das Problem, dass sich der
R
Code vom Python-Code unterscheidet, ist jedoch immer noch ein Rätsel ...quelle
R
undPython
etwas mehr untersuchen.Ich werde den R-Code kommentieren:
Sie setzen Variablen an den falschen Stellen zurück, dh die Variablen
min_mse
solltenInf
außerhalb derfor
Schleifeoptimal_alpha
initialisiert undNULL
dort initialisiert werden. Dies wird:Die Ausgabe sollte nun konsistent die kleinsten Alpha-Werte sein, da die Prädiktoren eine starke Kolinearität aufweisen und die Antwort nur aus den verfügbaren Prädiktoren besteht, dh es gibt keine redundante Variable, die LASSO auf Null setzen soll, in diesem Fall wir Ich möchte keine Regularisierung durchführen, dh das Kleinste
alpha
sollte das Beste sein. Sie können die Wirkung von MSE hier sehen:Beachten Sie, dass ich 50 Alphas auf derselben Skala wie Sie verwende. Bei einem Alpha-Index von 35 werden beide Variablen auf Null gesetzt, was bedeutet, dass das Modell immer dasselbe tut und die mse stagniert.
Ein besseres Problem, um MSE, CV und LASSO zu studieren
Das obige Problem ist für den LASSO nicht sehr interessant. Das LASSO führt die Modellauswahl durch, daher möchten wir sehen, dass es tatsächlich die interessierenden Parameter auswählt. Es ist beeindruckender zu sehen, dass das Modell tatsächlich ein Alpha auswählt, das die MSE tatsächlich senkt, dh uns bessere Vorhersagen gibt, indem einige Variablen weggeworfen werden. Hier ist ein besseres Beispiel, in dem ich eine Reihe redundanter Prädiktoren hinzufüge.
Dann führen Sie einfach die for-Schleife wie im obigen Code aus! Beachten Sie, dass ich das Maximum
alphas
von 6 von 6 auf 1,5 gesetzt habe, um den Effekt in der folgenden Darstellung zu sehen. Jetzt ist der bestealpha
Wert nicht der niedrigste, aber Sie können im Diagramm sehen, dass die Kreuzvalidierungs-MSE einen Tropfen nimmt und am Ende wieder ansteigt. Der niedrigste Punkt in diesem Diagramm entspricht dem Alpha-Index mit dem niedrigsten CV-Fehler.quelle