逻辑回归原理及实现

1.逻辑回归

逻辑回归是一种线性回归模型,它假设数据服从伯努力分布(二项分布,0-1分布),通过极大似然估计,运用梯度下降方法(牛顿法) 求解,进而达到二分类目的。逻辑回归与线性回归有很多相似之处，去除Sigmoid映射函数的话，逻辑回归算法就是一个线性回归。逻辑回归以线性回归理论作为支持。由于引入了Sigmoid函数，可以处理非线性问题，因此可以轻松处理0/1分布问题。

2.伯努利(Binomial)分布

伯努利分布又称为0-1分布、两点分布，通常对应事物的正反两面（成功与失败，是与否，正面与反面），指的是对于随机变量X有, 参数为p(0<p<1)，如果它分别以概率p和1-p取1和0为值。期望EX= p,方差DX=p(1-p)。逻辑回归中是对数据进行二分类，可以记为0,1。假设事件发生的概率为P(x=1)=P，则伯努利分布概率函数（公式）可以用下式表示：

　　其中上式的x只能取0或者1，1表示事件发生，0表示事件不发生。

3.Sigmoid函数与逻辑回归函数

　　Sigmoid函数公式如下:

 

　　这个公式有很多好的特性，g(z)函数值范围是(0,1)区间，当z取0时，值为0.5；远离0的地方，函数值很快接近0或者1。

 

　　逻辑回归函数形式如下：

　　上式也叫逻辑回归预测函数，预测在x,θ确定的情况下，y=1的概率，其中：

　　若预测y=0的概率，则：

 

　　为了更加方便的表示每个样本的概率，将y=0与y=1进行整合，得到下式：

4.似然函数

　　统计学中，似然函数是一种关于统计模型参数的函数。给定输出x时，关于参数θ的似然函数L(θ|x)（在数值上）等于给定参数θ后变量X的概率：L(θ|x)=P(X=x|θ)。简单理解就是，通过给定一支的样本来反推模型参数。也就是寻找参数θ使得样本发生的概率尽可能大。所以对于离散型随机变量而言，假设x的分布率为P(x|θ)，样本集D上有m个样本，则D上的似然函数为：

　　上式右边表示X1到Xm的概率的乘积，我们是寻找参数θ，让模型尽可能准确率高，也就是认为所有事件发生的概率尽可能大，所以就是求当似然函数取最大值时的因变量θ。也就是似然函数一般要构造成凸函数，可导才好求极大值。

5.对数似然函数

　　代价函数（Cost Function）是整个样本集的平均误差，对所有损失函数值的平均。在逻辑回归中，最常用的代价函数是交叉熵(Cross Entropy)函数，或者说是对数函数，他们本质是一样的。通过上文我们可以得到逻辑回归的似然函数为：

 

　　所以根据该似然函数理论上就能求出参数θ了。但是，乘积的式子很复杂，不便于求解。将似然函数转换为对数似然函数（以e为底）如下：

　　其中：

　　此时，采用数学的方法求极值，可求得最大值对应的参数θ。

6.对数似然损失函数与梯度求导

　　在机器学习中一般采用梯度下降法进行求解。上面式子求的是最大值，梯度下降求得是最小值，故添加一个负号,m个样本的平均损失函数转换为下面这个式子：

　　式子中有个m分之一，表示用了m个样本数据作为一组，表示求损失函数的平均值！式子中当真实值y=1，预测概率越接近1，损失值越小与经验一致。真实值y=0时，预测概率越接近0损失越小。下面是借鉴别人的求梯度式子推到。

　　最后的结果最后那个x变量中的i表示第i个样本，j表示第j个特征。括号内为预测减去真实值。

7.参数更新

　　利用梯度进行所要求的θ更新，如下式：

　　其中α为学习率。

8.召回率、精确率

　　精确率是针对我们预测结果而言的，它表示的是预测为正的样本中有多少是真正的正样本。那么预测为正就有两种可能了，一种就是把正类预测为正类(TP)，另一种就是把负类预测为正类(FP)，也就是：

　　召回率是针对我们原来的样本而言的，它表示的是样本中的正例有多少被预测正确了。那也有两种可能，一种是把原来的正类预测成正类(TP)，另一种就是把原来的正类预测为负类(FN)。

9.阈值指定与置信区间

　　预测时，如果不指定概率阈值，则根据sigmoid判断当概率大于等于0.5，判断为y=1。实际根据需要指定不同的阈值，比如概率大于0.9时，才认为预测准确。关于概率阈值具体研究，本文暂不做介绍。置信区间展现的是这个参数的真实值有一定概率落在测量结果的周围的程度。置信水平表示区间估计的把握程度，置信区间的跨度是置信水平的正函数，即要求的把握程度越大，势必得到一个较宽的置信区间，这就相应降低了估计的准确程度。比如，为了更严谨，更有把握预测癌症，完全可以把概率阈值设置为0.98，但是满足这样的区间一般会小了很多，把握是大了，但是能诊断出是否患癌症晚期的病人少了。

10.说明

　　逻辑回归损失函数并没有采用均方差函数，当模型完全预估错误时（y=1, p=0; 或y=0， p=1），损失是1。预估正确时，损失是0。错误值离正确值的“距离”相对较小，区分度不大。另外，上面的损失函数相对θ并非是凸函数，而是有很多极小值（local minimum）的函数。因此，很多凸优化的算法（如梯度下降）无法收敛到全局最优点。

11.代码实现

　　纯python代码实现:

logi.py

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn import preprocessing as pp
import numpy.random

def sigmoid(z):#参数z可以是一个实数，也可以是一个ndarray数据类型，多维矩阵都可以，比如向量[1,2,3]，结果会是[0.73,0.88,0.95]
    return 1/(1+np.exp(-z))

#定义预测函数model,预测的是概率，范围在(0,1);参数X应该如[[1,2,3],[4,5,6],...],theta如[1,2,3]
def model(X,theta):
    return sigmoid(np.dot(X,theta.T))#其结果如[[1],[2],[3],...]

#m表示多少个样本，也就是行数。y表示因变量真实值0或者1,x表示自变量，当前有两列。theta如[1,2,3],y_numpyarray如[1,0,1],x_array如[[1,2,3],[4,5,6],...]
def loss(m,theta,y_numpyarray,x_array):
    left=y_numpyarray*np.log(model(x_array,theta)+0.005)#结果如[[1],[2],[3],...]
    right=(1-y_numpyarray)*np.log(1-model(x_array,theta)+0.005)#结果如[[1],[2],[3]]
    return np.sum((left+right))/(-m)#结果如np.sum([1 2])，结果为一个实数

def loss_L2(m,theta,y_numpyarray,x_array,lamda):
    theta2=theta[1:]
    left=y_numpyarray*np.log(model(x_array,theta)+0.005)#结果如[[1],[2],[3],...]
    right=(1-y_numpyarray)*np.log(1-model(x_array,theta)+0.005)#结果如[[1],[2],[3]]
    l2=(lamda/m)*(sum(theta2*theta2))#L2惩罚系数
    return np.sum((left+right))/(-m)+l2#结果如np.sum([1 2])，结果为一个实数

#计算第i个样本的第j列（本例中是第1,2,3列）的梯度（导数）
def gradient(X,y,theta):#批量梯度下降，使用的是整个样本的平均损失函数,X为[[1,2,3],[4,5,6],...],y为[0,1,1,0,...],theta为[0.1,0.2,0.3]，都为ndarray类型
    grad=np.zeros(theta.shape)#grad如[0. 0. 0.]
    error=(model(X,theta)-y)#按公式计算m个样本，结果如[[1],[2],[3],...]
    for j in range(len(theta)):#计算第j个theta梯度
        term=np.multiply(error,X[:,j])#这里的X[:,j]如[[1],[2],[3]]，每个元素为第j列
        grad[j]=np.sum(term.ravel())/len(X)#ravel()将元素展平为[1,2,3,4...],sum求和后是一个数
    return grad

def gradient_L2(X,y,theta,lamda):#批量梯度下降，使用的是整个样本的平均损失函数,X为[[1,2,3],[4,5,6],...],y为[0,1,1,0,...],theta为[0.1,0.2,0.3]，都为ndarray类型
    theta2 = theta[1:]
    grad=np.zeros(theta.shape)#grad如[0. 0. 0.]
    error=(model(X,theta)-y)#按公式计算m个样本，结果如[[1],[2],[3],...]
    for j in range(len(theta)):#计算第j个theta梯度
        term=np.multiply(error,X[:,j])#这里的X[:,j]如[[1],[2],[3]]，每个元素为第j列
        grad[j]=np.sum(term.ravel())/len(X)+(lamda*(sum(theta2*theta2))/len(X))    #ravel()将元素展平为[1,2,3,4...],sum求和后是一个数
    return grad

#参数更新，更新theta
def updateParameter(alpha,gradtheta,theta):#alpha表示学习率，比如0.005,grad为梯度,[gradtheta1,gradtheta2,gradtheta3]，包括所有theta的梯度,theta是一个向量[theta1,theta2,theta3]。
    theta=theta-alpha*gradtheta
    return theta

def shuffleData(data):
    data=data.iloc[:,:].values.astype(float)
    # np.random.shuffle(data)
    cols=data.shape[1]
    X=data[:,0:cols-1]
    y=data[:,cols-1:]
    return X,y

#停止迭代三种方式
#第一种方式直接指定迭代次数，比如2000次，得反复调试观看损失图，才能确定次数多少合适。
#第二种方式看梯度，当梯度接近0时，表明损失函数最小了，此时停止迭代。
#第三种方式看损失函数的值，当前一次与后一次的损失插值接近0时停止。
#这里以第三种方式为例
def iterFlag(beforeLoss,lastLoss,minThreshold):#比如minThreshold=0.0000001，最后一次与倒数第二次损失值比这还小就停止迭代
    if(np.fabs(beforeLoss-lastLoss)<minThreshold):
        return True
    else:
        return False

#画图
def lossPlot(XRange,lossList):
    fig, ax = plt.subplots(figsize=(12, 4))
    ax.plot(XRange,lossList, 'r')
    ax.set_xlabel('iterations')
    ax.set_ylabel('loss')
    ax.set_title("loss-iterations".upper())
    plt.show()

def readcsv(path):
    df=pd.read_csv(path)
    return df

def insertOnes(df,col,name):
    return df.insert(col, name, 1)

# def scaleArray(numpyarray):

def logisticPredict(x_testArray,theta):#预测函数
    return model(x_testArray,theta)

#prob为阈值，默认取0.5，以0.5进行分类划分
def classifier(predictY,prob=0.5):
    result=np.ones((len(predictY)))
    result[predictY<=prob]=0
    return result

def Precision(Y,predictY):#精确率,针对预测的结果,1代表正类
    temp = predictY[Y == predictY]
    TP = len(temp[temp == 1])
    temp2 = predictY[Y != predictY]
    FP = len(temp2[temp2 == 1])
    if (TP + FP) == 0:
        return 0
    else:
        return TP / (TP + FP)
def Recall(Y,predictY):#召回率
    temp = predictY[Y == predictY]
    TP = len(temp[temp == 1])
    temp2 = predictY[Y != predictY]
    FN = len(temp2[temp2 == 0])
    print("FN:", FN)
    if (TP + FN)==0:
        return 0
    else:
        return TP / (TP + FN)

乳腺癌症预测：

import pandas as pd
import numpy as np
from sklearn import preprocessing as pp
from logi import insertOnes, shuffleData, loss, iterFlag, lossPlot, gradient, updateParameter, logisticPredict, \
    classifier, Precision, Recall

if __name__=='__main__':
    print("1.开始读取数据")
    # 构造列标签名字
    column = ['Sample code number', 'Clump Thickness', 'Uniformity of Cell Size', 'Uniformity of Cell Shape',
              'Marginal Adhesion', 'Single Epithelial Cell Size', 'Bare Nuclei', 'Bland Chromatin', 'Normal Nucleoli',
              'Mitoses', 'Class']
    cancerDataFrame =pd.read_csv('breast-cancer-wisconsin.data',names=column)#如果没有离线数据集，也可以使用sklearn库读取这个数据集datasets.load_breast_cancer
    print(cancerDataFrame.shape[1])
    col=cancerDataFrame.shape[1]
    # 缺失值进行处理，这里暂时不考虑这步，没有缺失值
    #数据处理
    print(cancerDataFrame.head())
    #去掉没用的列，本例中是sample code number
    cancerDataFrame.drop(['Sample code number'], axis=1, inplace=True)  # axis参数默认为0，为1表示删除列
    print(cancerDataFrame.head())
    # 缺失值数据处理，看数据比如41行有？
    cancerDataFrame=cancerDataFrame.replace(to_replace='?',value=np.nan)
    cancerDataFrame=cancerDataFrame.dropna()
    #数据删除后的行数
    rows=cancerDataFrame.shape[0]#行数
    #加入一列，全为1
    insertOnes(cancerDataFrame, 0, "ones")
    #打乱数据
    X,y=shuffleData(cancerDataFrame)#X,y都是ndarray类型
    y[y==2]=0
    y[y==4]=1
    print(y)
    x_numpyarray=X
    #对X进行标准化，据说标准化不会影响logistic回归，这里暂且写上，不需要的话可以注释这两行
    x_numpyarray = pp.scale(X)#这里使用sklearn自带的标准化函数
    x_numpyarray[:, 0] = 1#归一化后，第一列(theta常数,公式中的最后那个theta)重新设置为1
    # insertOnes(cancerDataFrame, 0, "ones")
    print(y)
    #数据训练测试集分割，这里是9:1
    x_trainArray=x_numpyarray[0:int(rows*0.8)]
    y_trainArray=y[0:int(rows*0.8)]
    x_testArray=x_numpyarray[int(rows*0.8):]
    y_testArray=y[int(rows*0.8):]
    #循环迭代
    beforeLoss=-1
    lastLoss=0
    # theta=np.zeros(x_numpyarray.shape[1])#初始设置全为0好像不行
    theta=np.ones(x_numpyarray.shape[1])
    # theta =np.array([0.4,12,0.8,9,0.4,2,1,1,4,3])
    minThreshold=1E-4
    print(rows)
    lossList=[]
    alpha=0.0001#学习率
    i=0
    while True:
        currentLoss=loss(rows,theta,y_trainArray,x_trainArray)
        i=i+1
        lossList.append(currentLoss)
        if i==500:
            break
        #调用梯度计算函数
        gradtheta=gradient(x_trainArray, y_trainArray, theta)
        #参数更新
        theta=updateParameter(alpha, gradtheta, theta)
    XRange=np.arange(0,len(lossList))
    lossPlot(XRange,lossList)
    #预测，预测的是概率，使用逻辑回归预测函数，通过概率进行二分类
    predictP=logisticPredict(x_testArray,theta)#预测恶性的概率
    result=np.dot(x_testArray, theta.T)
    print(theta)
    result=classifier(predictP,0.5)
    # 精确率，召回率
    Y=y_testArray.ravel()
    precision = Precision(Y, result)
    recall=Recall(Y, result)
    print("精准率",precision)
    print("召回率",recall)

12.结果分析

　　代码中随着迭代次数的增多，精准率与召回率及损失函数都在变化，但是当迭代次数太大时，损失函数可以达到最低值，但预测结果都不对，表明模型过拟合，繁化误差太大，此时模型系数theta出现其中一个权重很大，其他几乎为0（导致了过拟合）。可以调整迭代参数与学习率，使得精准率与召回率较高的水平。另外可以使用正则化惩罚系数可以有效避免过拟合，有效避免某一个参数很大的情况出现。过拟合问题往往源自过多的特征。

　　i==30，学习率=0.0001且标准化，无正则化时：

 

 

 

　　i==100，学习率=0.0001且标准化，无正则化时：

 

 

 

　　当i==1000时，学习率=0.0001且标准化，无正则化时

 

 

 

　　可以看到模型过拟合了，其中第一个权重特别大，且是我们插入的常数项偏置，参数太大表明模型过拟合，其导数太大，拟合曲线不够平滑。如果在不考虑迭代次数的情况下（迭代很多次），可以减少特征数量或者引入正则项惩罚系数。

13.正则化（regularization）

　　添加正则化项的作用是防止模型过拟合，让模型获得抗噪声的能力，这将提升模型对未知样本的预测性能。机器学习中常用的正则化有两种，L1正则化与L2正则化。

L1会趋向于产生少量的特征，而其他的特征都是0。最优的参数值很大概率出现在坐标轴上，这样就会导致某一维的权重为0 ，产生稀疏权重矩阵。L2会选择更多的特征，这些特征权重都会接近于0，而不为0。当目标损失函数与L1,L2范数第一次相交时，得到的是加入了惩罚项后的最优解（不为单个目标函数的最优解）。

　　损失函数:

　　根据惯例，上式中惩罚系数只针对特征而言，常数项偏置不用管。为了求导方便，可以选择分母乘以2。

　　梯度求导：

 

　　参数更新和非正则化参数更新差不多。

　　代码见上面的或者看我的百度网盘链接。

14.回归与变量的共线性问题

如果模型预测效果始终不好，可以考虑变量与变量之间是否存在多重共线性的问题，如果存在可以根据相关性剔除，保留大的一个，再做回归，这里只做简单介绍。

补充

　　可以直接调用python中sklearn库的函数实现逻辑回归，另外y分类标签据说不一定是0或1，任意两类都行，比如本数据的2或者4，大家可以试一下，而且损失函数此时可以为负值。可能是必须设置为0和1，分别代表良性与恶性。另外，关于决策边界的问题，这里暂不详细讲解。在癌症方面，召回率更重要，我们并不想癌症恶性预测为良性。本文的正则化似乎没有发挥用处，不知道原因，另外癌症数据集数据近700条，只训练十几轮，损失函数接近最小值。而具有100条数据的录取数据集1000轮迭代后，损失函数才最小。关于theta的设置，癌症数据集初始化全0，好像不太行，录取数据集好像可行。

参考资料:

https://blog.csdn.net/iqdutao/article/details/109478633

https://blog.csdn.net/qq_37018566/article/details/83241637

https://www.jianshu.com/p/54ed63a7f816

https://blog.csdn.net/u012344939/article/details/104335742

https://blog.csdn.net/u010159842/article/details/46426909

https://blog.csdn.net/kdongyi/article/details/83932945

数据及代码：

　　数据及本文代码百度网盘地址链接：https://pan.baidu.com/s/1DojPSLLn-ghewEjUU8Soxw   提取码：6t90 

 

本文任何错误与不足之处，恳请指正，欢迎交流！