邏輯回歸是分類當中極為常用的手段，它屬于機率型非線性回歸，分為二分類和多分類的回歸模型。對于二分類的logistic回歸，因變量y隻有“是”和“否”兩個取值，記為1和0。假設在自變量x1,x2,……,xp，作用下，y取“是”的機率是p，則取“否”的機率是1-p。下面将對最為常用的二分類logistic回歸模型的原理以及應用進行介紹。（不想看原理的可以直接調至後半部分，有代碼示範）
sigmoid函數
在logistic回歸的二分類問題中，要用到的函數就是sigmoid函數。sigmoid函數非常簡單，他的表達式是
（當z為0的時候，函數值為0.5；随着z的增大，函數值逼近于1；随着z的減小，函數值逼近于0.）
生成sigmoid函數圖的代碼：
import numpy
import math
import matplotlib.pyplot as plt
def sigmoid(x):
a = []
for item in x:
a.append(1.0/(1.0 + math.exp(-item)))
return a
x = numpy.arange(-10, 10, 0.1)
y = sigmoid(x)
plt.plot(x,y)
plt.yticks([0.0, 0.5, 1.0])
plt.axhline(y=0.5, ls='dotted', color='k')
plt.show()
sigmoid函數很适合做我們剛才提到的二分類的分類函數。假設輸入資料的特征是(x0, x1, x2, ..., xn)，我們在每個特征上乘以一個回歸系數 (w0, w1, w2, ... , wn)，然後累加得到sigmoid函數的輸入z：
那麼，輸出就是一個在0~1之間的值，我們把輸出大于0.5的資料分到1類，把輸出小于0.5的資料分到0類。這就是Logistic回歸的分類過程。
基于最優化方法的最佳回歸系數的确定
由上，可知logistic回歸的大緻流程如下，我們要做的就是确定出最佳的w=(w0, w1, w2, ... , wn)。
損失函數與極大似然函數
在logistic回歸裡面，用極大似然法來求解模型參數。關于似然函數的概念可以參考kevinGao的部落格
http://www.cnblogs.com/kevinGaoblog/archive/2012/03/29/2424346.html先定義似然函數(每個樣本都認為是獨立的)：
根據似然函數的概念，令似然函數最大的那個機率就是最合理的。我們想最大化似然函數，為友善計算，是以，我們取對數
可知，當權向量 w使l(w)最大的時候,w最合理。
梯度上升法求參數
梯度上升法的基本思想是：要找到某函數的最大值，最好的方法就是沿着該函數的梯度方向搜尋。如果函數為f，梯度記為D，a為步長，那麼梯度上升法的疊代公式為：w：w+a*Dwf(w)。該公式停止的條件是疊代次數達到某個指定值或者算法達到某個允許的誤差範圍。首先對對數的函數的梯度進行計算:
通過矩陣乘法直接表示成梯度:
設步長為α, 則疊代得到的新的權重參數為:
這樣我們通過梯度上升法做極大似然估計來做Logistic回歸的過程就很清楚了，剩下的我們就需要通過代碼來實作Logistic回歸吧。
代碼實作
資料集：學生的gre，gpa和rank資訊作為變量，預測是否admit，若admit=1代表錄取，admit=0代表不錄取。
import pandas as pd
import statsmodels.api as sm
import pylab as pl
import numpy as np
df = pd.read_csv("binary.csv")
# 浏覽資料集
print (df.head())
# admit gre gpa rank
#0 0 380 3.61 3
#1 1 660 3.67 3
#2 1 800 4.00 1
#3 1 640 3.19 4
#4 0 520 2.93 4
# 重命名'rank'列，因為dataframe中有個方法名也為'rank'
df.columns = ["admit", "gre", "gpa", "prestige"]
#資料統計情況
print (df.describe())
# admit gre gpa prestige
#count 400.000000 400.000000 400.000000 400.00000
#mean 0.317500 587.700000 3.389900 2.48500
#std 0.466087 115.516536 0.380567 0.94446
#min 0.000000 220.000000 2.260000 1.00000
#25% 0.000000 520.000000 3.130000 2.00000
#50% 0.000000 580.000000 3.395000 2.00000
#75% 1.000000 660.000000 3.670000 3.00000
#max 1.000000 800.000000 4.000000 4.00000
# 頻率表，表示prestige與admin的值相應的數量關系
print (pd.crosstab(df['admit'], df['prestige'], rownames=['admit']))
#prestige 1 2 3 4
#admit
#0 28 97 93 55
#1 33 54 28 12
拟變量（啞變量）
虛拟變量，也叫啞變量，可用來表示分類變量、非數量因素可能産生的影響。在計量經濟學模型，需要經常考慮屬性因素的影響。例如，職業、教育程度、季節等屬性因素往往很難直接度量它們的大小。隻能給出它們的“Yes—D=1”或”No—D=0”，或者它們的程度或等級。為了反映屬性因素和提高模型的精度，必須将屬性因素“量化”。通過構造0-1型的人工變量來量化屬性因素。pandas提供了一系列分類變量的控制。我們可以用get_dummies來将”prestige”一列虛拟化。
# 将prestige設為虛拟變量
dummy_ranks = pd.get_dummies(df['prestige'], prefix='prestige')
print (dummy_ranks.head())
# prestige_1 prestige_2 prestige_3 prestige_4
#0 0 0 1 0
#1 0 0 1 0
#2 1 0 0 0
#3 0 0 0 1
#4 0 0 0 1
建構需要進行邏輯回歸的資料框：
# 除admit、gre、gpa外，加入了上面常見的虛拟變量（注意，引入的虛拟變量列數應為虛拟變量總列數減1，減去的1列作為基準）
cols_to_keep = ['admit', 'gre', 'gpa']
data = df[cols_to_keep].join(dummy_ranks.ix[:, 'prestige_2':])
print (data.head())
# admit gre gpa prestige_2 prestige_3 prestige_4
#0 0 380 3.61 0 1 0
#1 1 660 3.67 0 1 0
#2 1 800 4.00 0 0 0
#3 1 640 3.19 0 0 1
#4 0 520 2.93 0 0 1
# 需要自行添加邏輯回歸所需的intercept變量
data['intercept'] = 1.0
根據上述的資料框執行邏輯回歸：
# 指定作為訓練變量的列，不含目标列`admit`
train_cols = data[data.columns[1:]]
# sigmoid函數
def sigmoid(inX): #sigmoid函數
return 1.0/(1+np.exp(-inX))
#梯度上升求最優參數
def gradAscent(dataMat, labelMat):
dataMatrix=np.mat(dataMat) #将讀取的資料轉換為矩陣
classLabels=np.mat(labelMat).transpose() #将讀取的資料轉換為矩陣
m,n = np.shape(dataMatrix)
alpha = 0.00001 #設定梯度的閥值，該值越大梯度上升幅度越大
maxCycles = 300 #設定疊代的次數，一般看實際資料進行設定，有些可能200次就夠了
weights = np.ones((n,1)) #設定初始的參數，并都賦預設值為1。注意這裡權重以矩陣形式表示三個參數。
for k in range(maxCycles):
h = sigmoid(dataMatrix*weights)
error = (classLabels - h) #求導後內插補點
weights = weights + alpha * dataMatrix.transpose()* error #疊代更新權重
return weights
#得到權重
weights=gradAscent(train_cols, data['admit']).getA()
#print (weights)
根據拟合出來的模型，可以進行預測：
# 在這邊為友善，我們将訓練集拷貝一份作為預測集（不包括 admin 列）
import copy
test_data = copy.deepcopy(data)
# 預測集也要添加intercept變量
test_data['intercept'] = 1.0
# 資料中的列要跟預測時用到的列一緻
predict_cols = test_data[test_data.columns[1:]]
# 進行預測，并将預測評分存入 predict 列中
predict=[]
test=np.mat(predict_cols)
for i in test:
sum=sigmoid(i*np.mat(weights))
print (sum)
if sum <= 0.5:
predict.append('0')
else:
predict.append('1')
test_data['predict']=predict
#計算預測準确率
predict_right=0
for i in range(0,400):
if int(test_data.loc[i,'admit'])==int(test_data.loc[i,'predict']):
predict_right=1+predict_right
else:
predict_right=predict_right
print ("預測準确率:")
print ("%.5f" %(predict_right/400))
#預測準确率:
#0.68250
由上，可知模型預測的準确率為68.25%，但往往我們會改進梯度上升方法以提高預測準确率，比如，改為随機梯度上升法。随機梯度上升法的思想是，每次隻使用一個資料樣本點來更新回歸系數。這樣就大大減小計算開銷。
def stocGradAscent(dataMatrix,classLabels):
m,n=shape(dataMatrix)
alpha=0.01
weights=ones(n)
for i in range(m):
h=sigmoid(sum(dataMatrix[i] * weights))#數值計算
error = classLabels[i]-h
weights=weights + alpha * error * dataMatrix[i] #array 和list矩陣乘法不一樣
return weights
同時，還可以改進随機梯度上升法，如下：
def stocGradAscent1(dataMatrix,classLabels,numIter=150):
m,n=shape(dataMatrix)
weights=ones(n)
for j in range(numIter):
dataIndex=list(range(m))
for i in range(m):
alpha=4/(1+i+j)+0.01#保證多次疊代後新資料仍然具有一定影響力
randIndex=int(random.uniform(0,len(dataIndex)))#減少周期波動
h=sigmoid(sum(dataMatrix[randIndex] * weights))
error=classLabels[randIndex]-h
weights=weights + alpha*dataMatrix[randIndex]*error
del(dataIndex[randIndex])
return weights
結語
由上，大家一定對logistic回歸有了一定的了解，如果你不願意自己定義函數調整參數，你也可以調用已有的包來進行logistic回歸分類，比如sklearn庫中的LogisticRegression，以及statsmodels庫中的Logit。
原文釋出時間為：2018-07-26
本文作者：胡蘿蔔醬
本文來自雲栖社群合作夥伴“
Python愛好者社群”，了解相關資訊可以關注“
”