python中的scikit-learn库提供了一系列机器学习的算法,其中包括Logistic Regression。不过默认的,该库提供的是加入L2惩罚项的Logistic回归。但是程序设置上,python与R语言中有点出入,使得默认设置下同样数据出现的结果相差很大。

R语言

在R语言中执行带惩罚项的Logistic回归,需要加载glmnet包。包中Logistic回归的优化目标是:

其中$\text{nobs}$代表样本量。惩罚项如果是L2惩罚,则是$\frac{1}{2} \parallel \theta \parallel_2^2$,其中不包含截距项。如果是L1惩罚,就是$ \theta _1$.

glmnet的设置:

  • family=‘binomial’ 代表Logistic回归;
  • alpha 表示elasticnet的混合参数:

所以$\alpha=1$表示lasso的L1惩罚,$\alpha=0$表示岭回归的L2惩罚; lambda 就是优化目标中的参数λ 默认的,glmnet会添加截距项,并对自变量进行标准化。

python

python的scikit-learn中,Logistic回归是参考的这篇文献。文献中二阶惩罚优化的目标函数是:

与R语言的区别在于,这里的惩罚项里包含截距项。另外参数设置上,这里的$C$相当于$\lambda$,都是作为一个权衡似然函数和惩罚项的参数。

若想要与R语言的优化目标一致,则需要满足:

  • 不含截距项
  • $n\lambda = \frac{1}{C}$

其中,$n$表示样本量,$\lambda$是R语言glmnet函数中的lambda参数,$C$是python中的参数。

Code

选择iris数据的前100行,作一个二分类的Logistic回归。设置python中$C=1$,R语言中应该有$n\lambda=1$,所以lambda应该是$\frac{1}{n}=0.01$。

无截距项

python:无截距项

from sklearn import linear_model
from sklearn import datasets
iris = datasets.load_iris()
index = iris.target != 2
logit = linear_model.LogisticRegression(C=1.0,  tol=1e-8, fit_intercept= False)
logit.fit(iris.data[index, ], iris.target[index])
print [logit.coef_, logit.intercept_]
[array([[-0.44234418, -1.48544863,  2.23987714,  1.01147778]]), 0.0]

R语言:无截距,不标准化

library(glmnet)
index = 1:100
iris.x = as.matrix(iris[index, 1:4])
iris.y = as.factor(as.character(iris[index, 5]))
# 不标准化自变量,不添加截距项
logit = glmnet(iris.x, iris.y, family="binomial", alpha=0, lambda=1/100, standardize=F, thresh =1e-8, intercept = F)
coef(logit)
5 x 1 sparse Matrix of class "dgCMatrix"
                     s0
(Intercept)   .        
Sepal.Length -0.4385104
Sepal.Width  -1.4868544
Petal.Length  2.2384771
Petal.Width   1.0031797

结果基本一致。

添加截距项?

最后,若想添加截距项后两个程序结果一致,可以修改python中的intercept_scaling参数,其参数解释为:

when self.fit_intercept is True, instance vector x becomes [x, self.intercept_scaling], i.e. a “synthetic” feature with constant value equals to intercept_scaling is appended to the instance vector. The intercept becomes intercept_scaling * synthetic feature weight Note! the synthetic feature weight is subject to l1/l2 regularization as all other features. To lessen the effect of regularization on synthetic feature weight (and therefore on the intercept) intercept_scaling has to be increased.

具体的计算原理不清楚,但依据解释,只要增大intercept_scaling就可以减少惩罚项对截距的影响。

先来看看python中添加截距项的默认结果:

logit = linear_model.LogisticRegression(C=1.0,  tol=1e-8, fit_intercept= True)
logit.fit(iris.data[index, ], iris.target[index])
print [logit.coef_, logit.intercept_]
[array([[-0.4070443 , -1.46126494,  2.23984278,  1.00849909]]), array([-0.26042082])]

增大intercept_scaling=100000的结果:

logit = linear_model.LogisticRegression(C=1.0,  tol=1e-8, fit_intercept= True, intercept_scaling=1e5)
logit.fit(iris.data[index, ], iris.target[index])
print [logit.coef_, logit.intercept_]
[array([[-0.44234418, -1.48544863,  2.23987714,  1.01147778]]), 0.0]

R语言中添加截距项的结果:

logit = glmnet(iris.x, iris.y, family="binomial", alpha=0, lambda=1/100, 
               standardize=F, thresh = 1e-8, intercept = T)
coef(logit)
5 x 1 sparse Matrix of class "dgCMatrix"
                     s0
(Intercept)  -6.6114025
Sepal.Length  0.4403477
Sepal.Width  -0.9070024
Petal.Length  2.3084749
Petal.Width   0.9623247

增大intercept_scaling=100000后与R语言的结果基本一致了。