逻辑回归：原理+代码

（作者：陈玓玏）

逻辑回归算是传统机器学习中最简单的模型了，它的基础是线性回归，为了弄明白逻辑回归，我们先来看线性回归。

一、线性回归

假设共N个样本，每个样本有M个特征，这样就产生了一个N*M大小的样本矩阵。令矩阵为X，第i个样本为Xi，第i个样本的第j个特征为Xij。令样本的观测向量为Y，第i个样本的观测值为Yi，那么就会有以下公式：
（X+[1]N*1）*W = Y
也就是说，对于一批已经存在的样本，我们通过其X和Y之间的关系，求得一个系数矩阵W，对于那些Y值未知的变量，我们通过X和W求得其可能的Y，也就是对其进行预测。

线性回归的合理性是已经证明过的，那就是任何变量之间的关系都是可以通过多项式进行拟合的，线性回归其实就是在求解一个多元一次方程。

二、逻辑回归

逻辑回归和线性回归很像，差别在于逻辑回归在线性回归的基础上加了一个激活函数。逻辑回归所使用的激活函数也分两种，一个是sigmoid函数，一个是softmax函数，前者用于二分类，称为逻辑回归，后者用于多分类，一般称为逻辑分类。线性回归输出的y值通过sigmoid函数，会产生N*1维的大小均在0和1之间的向量，也就是每个样本属于正样本的概率值，而y值通过softmax函数，会输出N*C维的矩阵，C表示样本预测值可能的类别数。

为什么要在线性回归的基础上用这两个激活函数呢？我的理解是，其一，算是对y值进行了一个归一化，因为单纯的线性回归，如果不对x进行归一化，那么拟合出来的y值可能会差别很大，那么在分类问题上，你阈值取多少才合适呢？这个范围就很大了，而sigmoid函数和softmax函数能够把y值映射到[0,1]区间上，表示样本属于正样本的概率大小，这个问题就解决了。其二，这两个函数都具有近似阶跃函数的特性，也就是能够把样本区分得足够开，这样能够让分类的准确率更高些，否则y的线性回归预测值靠近0时，很容易错分。

三、逻辑回归应用场景

逻辑回归适用于：
1、变量少的场景；
2、要求模型解释性强的场景。

适用于变量少的场景是因为逻辑回归基于线性回归，那么在用梯度下降法求解系数矩阵时，要求x矩阵是满秩的，也就是行列式不为0，这样才能求得唯一解，而变量很多的情况下，变量间很容易产生强相关，如果你使用sklearn，相关变量很多的情况下就会报singular error的错误。

适用于模型解释解释性强的场景是因为逻辑回归的本质其实是给每个变量分配了一个权重，那么当你模型结果产生错误时，你能够很容易地分析出是哪个变量导致了错误，而决策树这类的算法相对而言就要难分析一些，尤其是集成学习算法，因为它的建模过程、迭代过程没有那么明朗，也没有那么容易分析。

鉴于第二点因素，逻辑回归常常用于风控建模，因为风控建模对变量和模型的解释性要求极高。

四、逻辑回归代码

下面这段代码跑出来结果的正确率是0.8，还是可以接受的，嘿嘿。

import numpy as np
import pandas as pd
import math


#初始化权重矩阵
def initWeight(data):
    # x = len(data)
    y = len(data[0])
    w = np.ones((y,1))*1
    w = np.mat(w)
    print(w)

    return w

def gradOptimize(data_x,data_y,w,minLoss,stepLen,maxTurn,turn):
    turn = turn+1
    loss = 0
    y_val = []
    y_val = np.dot(data_x,w)
    # stepLen = stepLen / 1.1
    # print(y_val)
    # for j in range(len(y_seperate)):
    #     for i in range(len(y_seperate[0])):
    #         y_val.append(y_seperate[j,i]+y_seperate[j,i+1])
    # print(y_val.shape)
    for i in range(y_val.shape[0]):
        #print(i)
        p = 1/(1+math.exp(-y_val[i]))
        y_val[i] = p
        # loss = loss+i*math.log(p,2)+(1-i)*math.log((1-p),2)
        loss = loss + i * math.log(p, 2) + (1 - i) * math.log((1 - p), 2)
    # print(y_val.shape[0])
    loss = -loss/y_val.shape[0]
    # print(y_val)
    # print(loss)
    # print(w.shape)
    if abs(loss)>minLoss and turn<maxTurn:
        w = w-stepLen*np.dot(data_x.T,y_val-data_y)
        # print(w)
        # for i in range(w.shape[0]):
        #     for j in range(w.shape[1]):
        #         stepLen = stepLen / 1.1
        #         w[i,j] = w[i,j]-stepLen*(1/(1+math.exp(-i)))*(1-1/(1+math.exp(-i)))
        w, y_val = gradOptimize(data_x,data_y,w,minLoss,stepLen,maxTurn,turn)
        return w, y_val
    else:
        print('===============')
        print(y_val)
        return w,y_val

def oneScale(data):
    newdata = pd.DataFrame(data)
    for i in range(newdata.shape[1]-1):
        newdata.ix[:,i] = newdata.ix[:,i].apply(lambda x:x/abs(newdata.ix[:,i]).max())
    data = newdata.values
    return data

if __name__=='__main__':

    data=[
        [-0.017612,14.053064,0],
        [-0.752157,6.538620,0],
        [-1.322371,7.152853,0],
        [0.423363,11.054677,0],
        [0.569411,9.548755,0],
        [-0.026632,10.427743,0],
        [0.667394,12.741452,0],
        [1.347183,13.175500,0],
        [1.217916,9.597015,0],
        [-0.733928,9.098687,0],
        [1.416614,9.619232,0],
        [1.388610,9.341997,0],
        [0.317029,14.739025,0],
        [-0.576525,11.778922,0],
        [-1.781871,9.097953,0],
        [-1.395634,4.662541,1],
        [0.406704,7.067335,1],
        [-2.460150,6.866805,1],
        [0.850433,6.920334,1],
        [1.176813,3.167020,1],
        [-0.566606,5.749003,1],
        [0.931635,1.589505,1],
        [-0.024205,6.151823,1],
        [-0.036453,2.690988,1],
        [-0.196949,0.444165,1],
        [1.014459,5.754399,1],
        [1.985298,3.230619,1],
        [-1.693453,-0.557540,1],
        [-0.346811,-1.678730,1],
        [-2.124484,2.672471,1]
    ]
    #对变量值做归一化
    data = oneScale(data)
    X = np.mat(data.copy())[:, 0:-1]
    data_y = np.mat(data.copy())[:, -1]
    b = np.ones(data_y.shape)  # 添加全1列向量代表b偏量
    data_x = np.column_stack((b, X))  # 特征属性集和b偏量组成x
    # data_x = oneScale(data_x)
    # print(data_x)
    #print(data_x)

    # print(data_x,data_y)
    w = initWeight(data)
    w,y_val = gradOptimize(data_x,data_y,w,0.1,0.02,300,0)
    print(w)
    print(y_val)
    #print(np.multiply(data_x,w)[-1])

    # data = pd.DataFrame(np.dot(data_x,w),columns=['predict_y'])
    # data['predict_y'] = data['predict_y'].apply(lambda x:1/(1+math.exp(-x)))
    data = pd.DataFrame(y_val,columns=['predict_y'])
    # data['predict_y'] = y_val
    data['real_y'] = data_y
    correct_ratio = data[abs(data['predict_y']-data['real_y'])<=0.5].count()/data.count()
    print(data)
    print(correct_ratio)