逻辑回归

3P-2L

梯度~反向传播~实例~局部~凸函数

一、逻辑回归简介

　一句话概括：逻辑回归假设数据服从伯努利分布，通过极大化似然函数的方法，运用梯度下降来求解参数，来达到将数据二分类的目的

　这里面其实包含了5个点 1：逻辑回归的假设，2：逻辑回归的损失函数，3：逻辑回归的求解方法，4：逻辑回归的目的，5:逻辑回归如何分类。这些问题是考核你对逻辑回归的基本了解

• 逻辑回归的假设：

　　　逻辑回归的第一个基本假设是假设数据服从伯努利分布。伯努利分布有一个简单的例子是抛硬币，抛中为正面的概率是p，抛中为负面的概率是1-p，逻辑回归需要将线性模型进行一下映射，从而能用于分类

　　　这里的映射函数或者叫分类器叫做sigmoid函数，通过sigmoid函数分类器，逻辑回归的第二个假设是假设样本为正的概率是：

( P(y=1|x; heta ) = frac{1}{1+e^{-z}} = frac{1}{1+e^{- heta^Tx}} )

• 逻辑回归的损失函数：

　　　常用的损失函数有以下几种：

• • 0-1损失函数：

(L(Y,f(X)) = left{egin{matrix}1 & Y eq f(X)\ 0 & Y = f(X)end{matrix} ight.)

 

• • 平方损失函数：

(L(Y,f(X)) = (Y - f(X))^2)

 

• • 绝对损失函数：

(L(Y,f(X)) = | Y - f(X))|)

 

• • 对数损失函数或对数似然损失函数：

(L(Y,P(Y|X)) = -logP(Y|X))

 

　　　对于逻辑回归模型，使用的是对数损失函数作为代价函数，则本例中，逻辑回归的损失函数为：

(cost(y,p(y|x)) = left{egin{matrix}-logp(y|x) & if(y = 1)\ -log(1-p(y|x))& if(y = 0)end{matrix} ight.)

 

　　　将上面的两个表达式合并，则得到单个数据点上的log损失为：

(cost(y,p(y|x)) = -ylogp(y|x) - (1-y)log(1-p(y|x)))

　　　因为y只有两种取值情况，1或0，分别令y=1或y=0，即可得到原来的分段表达式，即两者是等价的。全体样本的损失函数则可表达为：

(cost(y,p(y|x)) = sum_{i=1}^{m} -y_ilogp(y_i|x_i) - (1-y_i)log(1-p(y_i|x_i)))

 

　　　其中(p(y|x))由前面的逻辑回归模型定义，令：

(p(y|x) = h_ heta(x) = frac{1}{1+e^{- heta^Tx}})

 

　　　

　　　则最终的损失函数为：

(cost(y,p(y|x)) = sum_{i=1}^{m} (-y_ilogfrac{1}{1+e^{- heta^Tx}} - (1-y_i)log(1-frac{1}{1+e^{- heta^Tx}})))

　　　即逻辑回归的损失函数是它的极大似然函数：

　

•  逻辑回归的求解方法：

　　　由于该极大似然函数无法直接求解，一般通过对该函数进行梯度下降来不断逼近最优解。因为就梯度下降本身来看的话就有随机梯度下降，批梯度下降，small batch 梯度下降三种方式

　　

• 逻辑回归的目的：

　　　将数据二分类，提高准确率

• 逻辑回归如何分类：

　　　逻辑回归作为回归，如何应用到分类上去？y值确实是一个连续的变量。逻辑回归的做法是划定一个阈值，y值大于这个阈值的是一类，y值小于这个阈值的是另外一类。阈值具体如何调整根据实际情况选择。一般会选择0.5做为阈值来划分。

二、Sigmoid函数

　如果我们使用线性回归来预测y的取值。这样做会导致线性回归的结果取值并不为0或1

　逻辑回归使用一个函数来归一化y值，使y的取值在区间(0,1)内，这个函数称为Logistic函数(logistic function)，也称为Sigmoid函数(sigmoid function)。函数公式如下：

( g(z) = frac{1}{1+e^{-z}} )

　函数图像如下：

　该函数导数有一个特性：

 

　经过Sigmoid函数转换后的y值有特殊的含义，它表示结果取1的概率，因此对于输入x分类结果为类别1和类别0的概率分别为：

( P(y=1|x; heta ) = h_ heta (x) )

( P(y=0|x; heta ) = 1 - h_ heta (x) )

 　对上面的表达式合并一下就是：

( P(y|x; heta ) = (h_ heta (x))^y(1-h_ heta (x))^{1-y} )

 

三、推导过程

　梯度上升

　　构建完似然函数，我们最大似然估计，最终推导出θ的迭代更新表达式

　　只不过这里用的不是梯度下降，而是梯度上升，因为这里是最大化似然函数不是最小化似然函数。我们假设训练样本相互独立，那么似然函数表达式为：

　　同样对似然函数取log，转换为：

 

　　转换后的似然函数对θ求偏导，在这里我们以只有一个训练样本的情况为例：

 

　　这样我们就得到了梯度上升每次迭代的更新方向，那么(  heta  )的迭代表达式为：

 

四、逻辑回归参数

　逻辑回归的模型调优一般考虑对其正则惩罚项系数进行调整，sigmoid函数的选择，在实际场景中一般是保证一定精确率（Precision）的情况下尽可能提高召回率（Recall rate）

　1）class_weight：类型权重参数，默认None，可以输入"balanced"或者一个字典

　　该参数用于标示分类模型中各种类型的权重，可以不输入，即不考虑权重，或者说所有类型的权重一样

　　如果选择输入的话，可以选择balanced让类库自己计算类型权重，或者我们自己输入各个类型的权重，比如对于0，1的二元模型，可以定义class_weight={ 0:0.9, 1:0.1 }，这样类型0的权重为90%，而类型1的权重为10%

　　如果class_weight选择balanced，那么类库会根据训练样本量来计算权重。某种类型样本量越多，则权重越低，样本量越少，则权重越高

　　sklearn的官方文档中，当class_weight为balanced时，类权重计算方法如下：

　　n_samples / (n_classes * np.bincount(y))，n_samples为样本数，n_classes为类别数量，np.bincount(y)会输出每个类的样本数，例如y=[1,0,0,1,1]，则np.bincount(y)=[2,3]

　　那么class_weight有什么作用呢？在分类模型中，我们经常会遇到两类问题：

　　第一种是误分类的代价很高。比如对合法用户和非法用户进行分类，将非法用户分类为合法用户的代价很高，我们宁愿将合法用户分类为非法用户，这时可以人工再甄别，但是却不愿将非法用户分类为合法用户。这时，我们可以适当提高非法用户的权重

　　第二种是样本是高度失衡的，比如我们有合法用户和非法用户的二元样本数据10000条，里面合法用户有9995条，非法用户只有5条，如果我们不考虑权重，则我们可以将所有的测试集都预测为合法用户，这样预测准确率理论上有99.95%

　　但是却没有任何意义。这时，我们可以选择balanced，让类库自动提高非法用户样本的权重

　　提高了某种分类的权重，相比不考虑权重，会有更多的样本分类划分到高权重的类别，从而可以解决上面两类问题

　2）sample_weight：样本权重参数

　　由于样本不平衡导致我们的模型预测能力下降。遇到这种情况，我们可以通过调节样本权重来尝试解决这个问题

　　调节样本权重的方法有两种，第一种是在class_weight使用balanced。第二种是在调用fit函数时，通过sample_weight来自己调节每个样本权重

　　在scikit-learn做逻辑回归时，如果上面两种方法都用到了，那么样本的真正权重是class_weight * sample_weight

　3）penalty：正则化参数：’l1’,’l2’.默认是’l2’

　　在调参时如果我们主要的目的只是为了解决过拟合，一般penalty选择L2正则化就够了。但是如果选择L2正则化发现还是过拟合，即预测效果差的时候，就可以考虑L1正则化

　　另外，如果模型的特征非常多，我们希望一些不重要的特征系数归零，从而让模型系数稀疏化的话，也可以使用L1正则化

　　penalty参数的选择会影响我们损失函数优化算法的选择。即参数solver的选择

　4）C：smaller values specify stronger regularization

五、手撕逻辑回归

　1）数据获取加header

import pandas as pd
import numpy as np
 
#数据没有标题，因此加上参数header
df = pd.read_csv('https://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/breast-cancer-wisconsin.data',
                 header=None)
 
column_names = ['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']
df.columns = column_names

df[df['Sample code number'] == 1057067]

　　我们发现有些特征存在空值：

　　去掉空值：

df = df.dropna(how='any')
print(df.shape) # 剩下683个样本：(683, 11)
df[df['Sample code number'] == 1057067]

　　可以发现已经不存在这个样本了，已经被去掉了：

　　分割数据：

import warnings
warnings.filterwarnings('ignore')
 
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

#一般1代表恶性，0代表良性（本数据集4恶性，所以将4变成1，将2变成0）
#df['Class'][data['Class'] == 4] = 1
#df['Class'][data['Class'] == 2] = 0
df.loc[ df['Class'] == 4, 'Class' ] = 1
df.loc[ df['Class'] == 2, 'Class' ] = 0

#Sample code number特征对分类没有作用，将数据集75%作为训练集，25%作为测试集
X_train, X_test, y_train, y_test = train_test_split(df[column_names[1:10]], df[column_names[10]],
                                                    test_size = 0.25, random_state = 33)

 
ss = StandardScaler()
X_train = ss.fit_transform(X_train)
X_test = ss.transform(X_test)

X_train = np.round(X_train, 2)
X_test = np.round(X_test, 2)

y_train = np.array(y_train)
y_test = np.array(y_test
                 )
print(X_train.shape, X_test.shape)  # (512, 9) ，(171, 9) 加和正好是683

　　展示一些数据吧：

print(X_train[0:5])
print('
')
print(X_test[0:5])
# 前五条
'''
[[-1.19 -0.69 -0.04 -0.62 -0.53 -0.67 -0.56 -0.59 -0.34]
 [-1.19  0.31 -0.04  2.55  0.37  1.87  0.69  1.09 -0.34]
 [ 0.24 -0.03 -0.72 -0.62 -0.53 -0.67 -0.97 -0.59 -0.34]
 [-1.19 -0.69 -0.72 -0.62 -0.53 -0.67 -0.14 -0.59 -0.34]
 [-1.19 -0.69 -0.72 -0.62 -0.53 -0.67 -0.14 -0.59 -0.34]]


[[ 0.24  0.31  0.65 -0.62  2.18 -0.67 -0.14  1.09 -0.34]
 [-0.47 -0.69 -0.72 -0.62 -0.53 -0.67 -0.56 -0.59 -0.34]
 [ 2.04  1.63  1.67 -0.27  2.18  1.87  0.27  1.77  4.97]
 [-0.12  1.63  0.99  0.44 -0.08  0.18  2.76  1.09 -0.34]
 [ 0.24 -0.36  0.3  -0.62 -0.98 -0.67 -0.97 -0.59 -0.34]]
'''

　　你要知道我们的目的就是找到w和b，能够让我们的训练集误差最小，而我们的数据集特征维度是：dim = 9，所以w长度就是9

　2）正向传播

　　我们先初始化一个w和b

def initialize_with_zeros(dim):
    w = np.zeros((dim,1))
    b = 0
    return w,b

# 这里dim = 9
'''
[[0.]
 [0.]
 [0.]
 [0.]
 [0.]
 [0.]
 [0.]
 [0.]
 [0.]]
'''

　　定义激励函数：

def sigmoid(z):
    a = 1/(1+np.exp(-z))
    return a

　　正向传播：

# 这里的X就是我们的训练集，只是形状为：(dim, 512)，也就是X_train进行了转置
# 因为机器学习领域，我们所指的向量都是列向量
# X_train.shape = (512, 9)，转置后为(9 , 512)，每一列就是一个样本的特征向量值

w, b = initialize_with_zeros(9)
X = X_train.T

A = sigmoid(np.dot(w.T,X)+b)
print(A.shape) # (1, 512)，得到这个向量，表示每个样本的wx + b的值经过激活函数
print(A)
# 因为初始化w的9个分量都是0，所以激活前值为0，经过激活函数得到0.5

其实这就是逻辑回归的第二个假设：假设该样本为正的概率：

( P(y=1|x; heta ) = frac{1}{1+e^{-z}} = frac{1}{1+e^{- heta^Tx}} )

 

　　损失：（该轮所有训练样本的平均损失：也就是上面的负对数损失函数）

m = X.shape[1] # 样本数目
Y = y_train
cost = -( np.sum( Y*np.log(A) + (1-Y)*np.log(1-A) ) ) / m
cost

'''
0.6931471805599453
'''

 　　有了该轮的损失，下面我们就是想法通过调整权重w和偏置b，重新计算cost，使得cost能够降低（比当前cost：0.69小）

 

　3）极大似然函数

　　上面我们也说过：逻辑回归的损失函数是它的极大似然函数（没毛病吧，老铁！我们只需要求一个最佳的参数w，能够使得所有样本的概率值之积最大）

　　可是这么多样本，概率都在0-1之间，相乘后的结果值会很小很小，所以才会取对数，求最大值（梯度下降来求最小值，所以我们可以对目标函数取负值，然后求新的目标函数最小值，就可以使用梯度下降啦）

　　我们可以参考上面的推导过程：我这里画了一个图，对与一个样本的梯度和所有样本的梯度计算以及如何更新一个权重$w_1$，到如何更新9个权重$w_9$，这里的512是训练样本数量，共9个特征值

　　自己计算可以一个样本一个样本来，程序可以利用矩阵运算一下子求出dw和db：

def propagate(w, b, X, Y):
    """
    传参:
    w -- 权重, shape： (dim, 1)
    b -- 偏置项, 一个标量
    X -- 数据集，shape： (dim, m),m为样本数
    Y -- 真实标签，shape： (1,m)

    返回值:
    cost， dw ，db，后两者放在一个字典grads里
    """
    #获取样本数m：
    m = X.shape[1]
    
    # 前向传播 ：
    A = sigmoid(np.dot(w.T,X)+b)    # 调用前面写的sigmoid函数    
    cost = -(np.sum(Y*np.log(A)+(1-Y)*np.log(1-A)))/m                 
    
    # 反向传播：
    dZ = A-Y
    dw = (np.dot(X,dZ.T)) / m
    db = (np.sum(dZ)) / m
  
    #返回值：
    grads = {"dw": dw,
             "db": db}
    
    return grads, cost

　　有了梯度，我们再根据学习率，来求出每次的步长，来进行梯度上升或者下降（上升或者下降取决于你对损失函数是否添加了符号）

　　然后我们经过多轮训练，反复调整权重$w$：

def optimize(w, b, X, Y, num_iterations, learning_rate, print_cost = False):
    #定义一个costs数组，存放每若干次迭代后的cost，从而可以画图看看cost的变化趋势：
    costs = []
    #进行迭代：
    for i in range(num_iterations):
        # 用propagate计算出每次迭代后的cost和梯度：
        grads, cost = propagate(w,b,X,Y)
        dw = grads["dw"]
        db = grads["db"]
        
        # 用上面得到的梯度来更新参数：
        w = w - learning_rate*dw
        b = b - learning_rate*db
        
        # 每100次迭代，保存一个cost看看：
        if i % 100 == 0:
            costs.append(cost)
        
        # 这个可以不在意，我们可以每100次把cost打印出来看看，从而随时掌握模型的进展：
        if print_cost and i % 100 == 0:
            print ("Cost after iteration %i: %f" %(i, cost))
            
    #迭代完毕，将最终的各个参数放进字典，并返回：
    params = {"w": w,
              "b": b}
    grads = {"dw": dw,
             "db": db}
    return params, grads, costs

　　最后我们就可以根据训练结果来进行预测啦：

def predict(w,b,X):
    m = X.shape[1]
    Y_prediction = np.zeros((1,m))

    A = sigmoid(np.dot(w.T,X)+b)
    
    for  i in range(m):
        if A[0,i]>0.5:
            Y_prediction[0,i] = 1
        else:
            Y_prediction[0,i] = 0

    return Y_prediction

def logistic_model(X_train,Y_train,X_test,Y_test,learning_rate=0.1,num_iterations=2000,print_cost=False):
    #获特征维度，初始化参数：
    dim = X_train.shape[1]
    W,b = initialize_with_zeros(dim)

    #梯度下降，迭代求出模型参数：
    params,grads,costs = optimize(W,b,X_train.T,Y_train,num_iterations,learning_rate,print_cost)
    W = params['w']
    b = params['b']

    #用学得的参数进行预测：
    prediction_train = predict(W,b,X_train.T)
    prediction_test = predict(W,b,X_test.T)
    
    #计算准确率，分别在训练集和测试集上：
    accuracy_train = 1 - np.mean(np.abs(prediction_train - Y_train))
    accuracy_test = 1 - np.mean(np.abs(prediction_test - Y_test))
    print("Accuracy on train set:",accuracy_train )
    print("Accuracy on test set:",accuracy_test )

   #为了便于分析和检查，我们把得到的所有参数、超参数都存进一个字典返回出来：
    d = {"costs": costs,
         "Y_prediction_test": prediction_test , 
         "Y_prediction_train" : prediction_train , 
         "w" : w, 
         "b" : b,
         "learning_rate" : learning_rate,
         "num_iterations": num_iterations,
         "train_acy":accuracy_train,
         "test_acy":accuracy_test
        }
    return d

d = logistic_model(X_train, np.array(y_train), X_test, np.array(y_test), learning_rate = 0.005, num_iterations = 2000,  print_cost = True)

　　这里矩阵运算注意来回转置，其实不必去注意这种转换，个人觉得没有必要去扣这种细节，反正原矩阵相乘不行，那就转置相乘，就到这里吧