使用 Iris 数据集在 Python 中再现 LASSO / Logistic 回归结果
Reproducing LASSO / Logistic Regression results in R with Python using the Iris Dataset
我正在尝试在 Python 中重现以下 R 结果。在这种特殊情况下，R 预测技能低于 Python 技能，但根据我的经验通常不是这种情况（因此想要重现 Python 中的结果），所以请忽略此处详细说明。
目的是预测花的种类（'versicolor' 0 or 'virginica' 1）。我们有 100 个标记样本，每个样本包含 4 个花特征：萼片长度、萼片宽度、花瓣长度、花瓣宽度。我将数据分成训练集（60% 的数据）和测试集（40% 的数据）。对训练集应用 10 折交叉验证来搜索最优的 lambda（优化的参数是 scikit-learn 中的"C"）。
我正在使用 glmnet in R with alpha set to 1 (for the LASSO penalty), and for python, scikit-learn's LogisticRegressionCV 函数和 "liblinear" 求解器（唯一可以与 L1 惩罚一起使用的求解器）。交叉验证中使用的评分指标在两种语言之间是相同的。然而，模型结果不知何故不同（为每个特征找到的截距和系数差异很大）。
R代码
library(glmnet)
library(datasets)
data(iris)
y <- as.numeric(iris[,5])
X <- iris[y!=1, 1:4]
y <- y[y!=1]-2
n_sample = NROW(X)
w = .6
X_train = X[0:(w * n_sample),] # (60, 4)
y_train = y[0:(w * n_sample)] # (60,)
X_test = X[((w * n_sample)+1):n_sample,] # (40, 4)
y_test = y[((w * n_sample)+1):n_sample] # (40,)
# set alpha=1 for LASSO and alpha=0 for ridge regression
# use class for logistic regression
set.seed(0)
model_lambda <- cv.glmnet(as.matrix(X_train), as.factor(y_train),
nfolds = 10, alpha=1, family="binomial", type.measure="class")
best_s <- model_lambda$lambda.1se
pred <- as.numeric(predict(model_lambda, newx=as.matrix(X_test), type="class" , s=best_s))
# best lambda
print(best_s)
# 0.04136537
# fraction correct
print(sum(y_test==pred)/NROW(pred))
# 0.75
# model coefficients
print(coef(model_lambda, s=best_s))
#(Intercept) -14.680479
#Sepal.Length 0
#Sepal.Width 0
#Petal.Length 1.181747
#Petal.Width 4.592025
Python代码
from sklearn import datasets
from sklearn.linear_model import LogisticRegressionCV
from sklearn.preprocessing import StandardScaler
import numpy as np
iris = datasets.load_iris()
X = iris.data
y = iris.target
X = X[y != 0] # four features. Disregard one of the 3 species.
y = y[y != 0]-1 # two species: 'versicolor' (0), 'virginica' (1). Disregard one of the 3 species.
n_sample = len(X)
w = .6
X_train = X[:int(w * n_sample)] # (60, 4)
y_train = y[:int(w * n_sample)] # (60,)
X_test = X[int(w * n_sample):] # (40, 4)
y_test = y[int(w * n_sample):] # (40,)
X_train_fit = StandardScaler().fit(X_train)
X_train_transformed = X_train_fit.transform(X_train)
clf = LogisticRegressionCV(n_jobs=2, penalty='l1', solver='liblinear', cv=10, scoring = ‘accuracy’, random_state=0)
clf.fit(X_train_transformed, y_train)
print clf.score(X_train_fit.transform(X_test), y_test) # score is 0.775
print clf.intercept_ #-1.83569557
print clf.coef_ # [ 0, 0, 0.65930981, 1.17808155] (sepal length, sepal width, petal length, petal width)
print clf.C_ # optimal lambda: 0.35938137
您在这里遇到的问题是数据集的排序（请注意，我没有检查 R 代码，但我确定这就是问题所在）。如果我 运行 你的代码然后 运行 这个
print np.bincount(y_train) # [50 10]
print np.bincount(y_test) # [ 0 40]
你可以看到训练集不代表测试集。但是，如果我对您的 Python 代码进行一些更改，那么我得到的测试准确度为 0.9.
from sklearn import datasets
from sklearn import preprocessing
from sklearn import model_selection
from sklearn.linear_model import LogisticRegressionCV
from sklearn.preprocessing import StandardScaler
import numpy as np
iris = datasets.load_iris()
X = iris.data
y = iris.target
X = X[y != 0] # four features. Disregard one of the 3 species.
y = y[y != 0]-1 # two species: 'versicolor' (0), 'virginica' (1). Disregard one of the 3 species.
X_train, X_test, y_train, y_test = model_selection.train_test_split(X, y,
test_size=0.4,
random_state=42,
stratify=y)
X_train_fit = StandardScaler().fit(X_train)
X_train_transformed = X_train_fit.transform(X_train)
clf = LogisticRegressionCV(n_jobs=2, penalty='l1', solver='liblinear', cv=10, scoring = 'accuracy', random_state=0)
clf.fit(X_train_transformed, y_train)
print clf.score(X_train_fit.transform(X_test), y_test) # score is 0.9
print clf.intercept_ #0.
print clf.coef_ # [ 0., 0. ,0., 0.30066888] (sepal length, sepal width, petal length, petal width)
print clf.C_ # [ 0.04641589]
我不得不对这里的几件事感到不满。
首先，"for python, scikit-learn's LogisticRegressionCV function with the "liblinear“求解器（唯一可以与L1惩罚一起使用的求解器）”。这显然是错误的，除非您打算以某种更明确的方式对其进行限定。只需看一下 sklearn.linear_model 类 的描述，您就会看到一些特别提到 L1 的内容。我相信其他人也允许您实施它，但我真的不想数数。
其次，您拆分数据的方法不太理想。看一下你拆分后的输入输出，你会发现在你的拆分中所有的测试样本的目标值都是1，而1的目标值只占你训练样本的1/6。这种不代表目标分布的不平衡将导致您的模型拟合不佳。例如，仅使用开箱即用的 sklearn.model_selection.train_test_split，然后完全按照原来的方式重新装配 LogisticRegressionCV 分类器，结果的准确度为 .92
综上所述，现在有一个 glmnet package for python，您可以使用此软件包复制您的结果。该项目的作者有一篇博客讨论了尝试使用 sklearn 重新创建 glmnet 结果的一些限制。具体来说：
"Scikit-Learn has a few solvers that are similar to glmnet, ElasticNetCV and LogisticRegressionCV, but they have some limitations. The first one only works for linear regression and the latter does not handle the elastic net penalty." - Bill Lattner GLMNET FOR PYTHON
上面的例子有几点不同：
系数的比例
- glmnet (https://cran.r-project.org/web/packages/glmnet/glmnet.pdf) 标准化数据和 "The coefficients are always returned on the original scale"。因此，您在调用 glmnet 之前没有缩放数据。
- Python 代码标准化数据，然后适合该标准化数据。本例中的系数采用标准化比例，而不是原始比例。这使得示例之间的系数不可比较。
LogisticRegressionCV 默认使用分层。 glmnet 使用 k-fold.
他们在拟合不同的方程。请注意，scikit-learn 逻辑适合 (http://scikit-learn.org/stable/modules/linear_model.html#logistic-regression) 与逻辑方面的正则化。 glmnet 将正则化放在惩罚上。
选择要尝试的正则化强度 - glmnet 默认为要尝试的 100 个 lambda。 scikit LogisticRegressionCV 默认为 10。由于 scikit 求解的方程，范围在 1e-4 和 1e4 (http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegressionCV.html#sklearn.linear_model.LogisticRegressionCV) 之间。
容忍度不同。在我遇到的一些问题中，收紧公差会显着改变系数。
- glmnet 默认 thresh 到 1e-7
- LogisticRegressionCV 默认 tol 到 1e-4
- 即使使它们相同，它们也可能测量的不是同一件事。我不知道 liblinear 措施是什么。 glmnet-"Each inner coordinate-descent loop continues until the maximum change in the objective after any coefficient update is less than thresh times the null deviance."
您可能想尝试打印正则化路径以查看它们是否非常相似，只是停在不同的强度上。然后你可以研究为什么。
即使更改了您可以更改的内容（并非以上所有内容），您也可能无法获得相同的系数或结果。尽管您在不同的软件中解决相同的问题，但软件解决问题的方式可能会有所不同。我们看到不同的尺度、不同的方程、不同的默认值、不同的求解器等。
我正在尝试在 Python 中重现以下 R 结果。在这种特殊情况下，R 预测技能低于 Python 技能，但根据我的经验通常不是这种情况（因此想要重现 Python 中的结果），所以请忽略此处详细说明。
目的是预测花的种类（'versicolor' 0 or 'virginica' 1）。我们有 100 个标记样本，每个样本包含 4 个花特征：萼片长度、萼片宽度、花瓣长度、花瓣宽度。我将数据分成训练集（60% 的数据）和测试集（40% 的数据）。对训练集应用 10 折交叉验证来搜索最优的 lambda（优化的参数是 scikit-learn 中的"C"）。
我正在使用 glmnet in R with alpha set to 1 (for the LASSO penalty), and for python, scikit-learn's LogisticRegressionCV 函数和 "liblinear" 求解器（唯一可以与 L1 惩罚一起使用的求解器）。交叉验证中使用的评分指标在两种语言之间是相同的。然而，模型结果不知何故不同（为每个特征找到的截距和系数差异很大）。
R代码
library(glmnet)
library(datasets)
data(iris)
y <- as.numeric(iris[,5])
X <- iris[y!=1, 1:4]
y <- y[y!=1]-2
n_sample = NROW(X)
w = .6
X_train = X[0:(w * n_sample),] # (60, 4)
y_train = y[0:(w * n_sample)] # (60,)
X_test = X[((w * n_sample)+1):n_sample,] # (40, 4)
y_test = y[((w * n_sample)+1):n_sample] # (40,)
# set alpha=1 for LASSO and alpha=0 for ridge regression
# use class for logistic regression
set.seed(0)
model_lambda <- cv.glmnet(as.matrix(X_train), as.factor(y_train),
nfolds = 10, alpha=1, family="binomial", type.measure="class")
best_s <- model_lambda$lambda.1se
pred <- as.numeric(predict(model_lambda, newx=as.matrix(X_test), type="class" , s=best_s))
# best lambda
print(best_s)
# 0.04136537
# fraction correct
print(sum(y_test==pred)/NROW(pred))
# 0.75
# model coefficients
print(coef(model_lambda, s=best_s))
#(Intercept) -14.680479
#Sepal.Length 0
#Sepal.Width 0
#Petal.Length 1.181747
#Petal.Width 4.592025
Python代码
from sklearn import datasets
from sklearn.linear_model import LogisticRegressionCV
from sklearn.preprocessing import StandardScaler
import numpy as np
iris = datasets.load_iris()
X = iris.data
y = iris.target
X = X[y != 0] # four features. Disregard one of the 3 species.
y = y[y != 0]-1 # two species: 'versicolor' (0), 'virginica' (1). Disregard one of the 3 species.
n_sample = len(X)
w = .6
X_train = X[:int(w * n_sample)] # (60, 4)
y_train = y[:int(w * n_sample)] # (60,)
X_test = X[int(w * n_sample):] # (40, 4)
y_test = y[int(w * n_sample):] # (40,)
X_train_fit = StandardScaler().fit(X_train)
X_train_transformed = X_train_fit.transform(X_train)
clf = LogisticRegressionCV(n_jobs=2, penalty='l1', solver='liblinear', cv=10, scoring = ‘accuracy’, random_state=0)
clf.fit(X_train_transformed, y_train)
print clf.score(X_train_fit.transform(X_test), y_test) # score is 0.775
print clf.intercept_ #-1.83569557
print clf.coef_ # [ 0, 0, 0.65930981, 1.17808155] (sepal length, sepal width, petal length, petal width)
print clf.C_ # optimal lambda: 0.35938137
您在这里遇到的问题是数据集的排序（请注意，我没有检查 R 代码，但我确定这就是问题所在）。如果我 运行 你的代码然后 运行 这个
print np.bincount(y_train) # [50 10]
print np.bincount(y_test) # [ 0 40]
你可以看到训练集不代表测试集。但是，如果我对您的 Python 代码进行一些更改，那么我得到的测试准确度为 0.9.
from sklearn import datasets
from sklearn import preprocessing
from sklearn import model_selection
from sklearn.linear_model import LogisticRegressionCV
from sklearn.preprocessing import StandardScaler
import numpy as np
iris = datasets.load_iris()
X = iris.data
y = iris.target
X = X[y != 0] # four features. Disregard one of the 3 species.
y = y[y != 0]-1 # two species: 'versicolor' (0), 'virginica' (1). Disregard one of the 3 species.
X_train, X_test, y_train, y_test = model_selection.train_test_split(X, y,
test_size=0.4,
random_state=42,
stratify=y)
X_train_fit = StandardScaler().fit(X_train)
X_train_transformed = X_train_fit.transform(X_train)
clf = LogisticRegressionCV(n_jobs=2, penalty='l1', solver='liblinear', cv=10, scoring = 'accuracy', random_state=0)
clf.fit(X_train_transformed, y_train)
print clf.score(X_train_fit.transform(X_test), y_test) # score is 0.9
print clf.intercept_ #0.
print clf.coef_ # [ 0., 0. ,0., 0.30066888] (sepal length, sepal width, petal length, petal width)
print clf.C_ # [ 0.04641589]
我不得不对这里的几件事感到不满。
首先，"for python, scikit-learn's LogisticRegressionCV function with the "liblinear“求解器（唯一可以与L1惩罚一起使用的求解器）”。这显然是错误的，除非您打算以某种更明确的方式对其进行限定。只需看一下 sklearn.linear_model 类 的描述，您就会看到一些特别提到 L1 的内容。我相信其他人也允许您实施它，但我真的不想数数。
其次，您拆分数据的方法不太理想。看一下你拆分后的输入输出，你会发现在你的拆分中所有的测试样本的目标值都是1，而1的目标值只占你训练样本的1/6。这种不代表目标分布的不平衡将导致您的模型拟合不佳。例如，仅使用开箱即用的 sklearn.model_selection.train_test_split，然后完全按照原来的方式重新装配 LogisticRegressionCV 分类器，结果的准确度为 .92
综上所述，现在有一个 glmnet package for python，您可以使用此软件包复制您的结果。该项目的作者有一篇博客讨论了尝试使用 sklearn 重新创建 glmnet 结果的一些限制。具体来说：
"Scikit-Learn has a few solvers that are similar to glmnet, ElasticNetCV and LogisticRegressionCV, but they have some limitations. The first one only works for linear regression and the latter does not handle the elastic net penalty." - Bill Lattner GLMNET FOR PYTHON
上面的例子有几点不同：
系数的比例
- glmnet (https://cran.r-project.org/web/packages/glmnet/glmnet.pdf) 标准化数据和 "The coefficients are always returned on the original scale"。因此，您在调用 glmnet 之前没有缩放数据。
- Python 代码标准化数据，然后适合该标准化数据。本例中的系数采用标准化比例，而不是原始比例。这使得示例之间的系数不可比较。
LogisticRegressionCV 默认使用分层。 glmnet 使用 k-fold.
他们在拟合不同的方程。请注意，scikit-learn 逻辑适合 (http://scikit-learn.org/stable/modules/linear_model.html#logistic-regression) 与逻辑方面的正则化。 glmnet 将正则化放在惩罚上。
选择要尝试的正则化强度 - glmnet 默认为要尝试的 100 个 lambda。 scikit LogisticRegressionCV 默认为 10。由于 scikit 求解的方程，范围在 1e-4 和 1e4 (http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegressionCV.html#sklearn.linear_model.LogisticRegressionCV) 之间。
容忍度不同。在我遇到的一些问题中，收紧公差会显着改变系数。
- glmnet 默认 thresh 到 1e-7
- LogisticRegressionCV 默认 tol 到 1e-4
- 即使使它们相同，它们也可能测量的不是同一件事。我不知道 liblinear 措施是什么。 glmnet-"Each inner coordinate-descent loop continues until the maximum change in the objective after any coefficient update is less than thresh times the null deviance."
您可能想尝试打印正则化路径以查看它们是否非常相似，只是停在不同的强度上。然后你可以研究为什么。
即使更改了您可以更改的内容（并非以上所有内容），您也可能无法获得相同的系数或结果。尽管您在不同的软件中解决相同的问题，但软件解决问题的方式可能会有所不同。我们看到不同的尺度、不同的方程、不同的默认值、不同的求解器等。