机器学习（三）：逻辑回归

一、问题的引入

在分类问题中，我们尝试预测的是结果是否属于某一个类(例如正确或错误)。比如:判断一封电子邮件是否是垃圾邮件;区别一个肿瘤是恶性的还是良性的。以肿瘤为例，用x代表各种检测指标，y代表肿瘤检验结果，为恶性则值为1，为良性则为0。这是一种离散输出，这就要引出逻辑回归算法了，该算法的输出值永远在0到1之间。

二、逻辑回归模型

1、函数：

注：x表示特征向量
g表示逻辑函数，是一个常用的逻辑函数为S形函数，公式为：将g的式子带入，于是：

2、根据图像，可知：

3、整合的式子


那接下来的问题就是如何拟合逻辑回归模型的参数?啦

3、代价函数

(1）引入似然函数，转化为对数似然，在线性回归要求解的是最小值，现在逻辑回归得到的是该事件发生的最大值。为了沿用梯度下降来求解，可以添加一个负号和常数解决。

（2）得到最终代价函数：

3、梯度下降求使得代价函数最小的参数

（1）求偏导

（2）得到结果

三、代码实现

1、前言

建立一个逻辑回归模型来预测一个学生是否被大学录取。有以前的申请人的历史数据，可以用它作为逻辑回归的训练集。对于每一个培训例子，你有两个考试的申请人的分数和录取决定。

2、数据

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import os
path = 'data' + os.sep + 'LogiReg_data.txt'
pdData = pd.read_csv(path, header=None, names=['Exam 1', 'Exam 2', 'Admitted'])
pdData.head()

得到结果：

Exam 1	        Exam 2	    Admitted
0	34.623660	78.024693	0
1	30.286711	43.894998	0
2	35.847409	72.902198	0
3	60.182599	86.308552	1
4	79.032736	75.344376	

看数据的维度

pdData.shape

得到结果

(100,3)

绘制散点图：

positive = pdData[pdData['Admitted'] == 1] # returns the subset of rows such Admitted = 1, i.e. the set of *positive* examples
negative = pdData[pdData['Admitted'] == 0] # returns the subset of rows such Admitted = 0, i.e. the set of *negative* examples

fig, ax = plt.subplots(figsize=(10,5))
ax.scatter(positive['Exam 1'], positive['Exam 2'], s=30, c='b', marker='o', label='Admitted')
ax.scatter(negative['Exam 1'], negative['Exam 2'], s=30, c='r', marker='x', label='Not Admitted')
ax.legend()
ax.set_xlabel('Exam 1 Score')
ax.set_ylabel('Exam 2 Score')

输出结果如下：

2、建立分类器

目标：建立分类器（求解三个参数），设定阀值，计算是否录取
(1)Sigmoid函数

import numpy as np
import pandas as pd 
import matplotlib.pyplot as plt
def sigmoid(z):
    return 1/(1+np.exp(-z))
# 建图
nums=np.arange(-10,10,step=1)
fig,ax=plt.subplots(figsize=(12,4))
ax.plot(nums,sigmoid(nums),"r")

（2）model函数
np.dot()是矩阵乘法，将数据传入sigmoid(z)中
np.dot(X,theta.T)相当于：

def model(X,theta):
    return sigmoid(np.dot(X,theta.T))

pdData.insert(0, 'Ones', 1)# 在第一列添加值全为1的列，列名是Ones
orig_data = pdData.as_matrix()#将frame转换为Numpy数组，
cols = orig_data.shape[1]
X = orig_data[:,0:cols-1]# 1到倒数第一列之前的数据
y = orig_data[:,cols-1:cols]# 倒数第一列的数据
theta = np.zeros([1,3])# 一行三列的矩阵全为0
print(X[:5])

输出结果如下：

[[  1.          34.62365962  78.02469282]
 [  1.          30.28671077  43.89499752]
 [  1.          35.84740877  72.90219803]
 [  1.          60.18259939  86.3085521 ]
 [  1.          79.03273605  75.34437644]]

print(y[0:5])

输出：

array([[ 0.],
       [ 0.],
       [ 0.],
       [ 1.],
       [ 1.]])

3、损失函数，计算损失

将对数似然函数去负号

$$D(h_{\theta }(x),y)=-y\log (h_{\theta }(x))-(1-y)\log (1-h_{\theta }(x))$$
求平均损失
$$J(\theta )=\frac{1}{n}\sum _{i=1}^{n}D(h_{\theta }(x_i),y_i)$$

import logging as log
def cost(X,y,theta):
    left = np.multiply(-y,np.log(model(X,theta)))
    right = np.multiply(1-y,np.log(1-model(X,theta)))
    return np.sum(left-right)/(len(X))
print(cost(X,y,theta))

得到结果:

0.69314718055994529

4、计算每个参数的梯度

$$\frac{\partial J}{\partial {\theta }_j}=-\frac{1}{m}\sum _{i=1}^n(y_i-h_{\theta }(x_i))x_{ij}$$

构造theta参数的时候，有三个参数，求解也应该有三个梯度

def gradient(X, y, theta):
    grad = np.zeros(theta.shape) # 初始化
    error = (model(X, theta)- y).ravel()#这里error和表达式反了，原因是把负号提到里面去了，revel()是将多维数组降为了一维
    for j in range(len(theta.ravel())): #分别对theta1、theta2做偏导
        term = np.multiply(error, X[:,j])
        grad[0, j] = np.sum(term) / len(X)
    
    return grad

三、三种不同的梯度下降方法

1、三种停止策略

（1）根据迭代次数进行停止。更新一次参数就是完成一次迭代
（2）根据迭代前后目标函数的变化进行停止，几乎没变化，就停下
（3）根据梯度，若前后迭代梯度的值差不多，就停下

STOP_ITER = 0
STOP_COST = 1
STOP_GRAD = 2

def stopCriterion(type, value, threshold):
    #设定三种不同的停止策略
    if type == STOP_ITER:        return value > threshold
    elif type == STOP_COST:      return abs(value[-1]-value[-2]) < threshold
    elif type == STOP_GRAD:      return np.linalg.norm(value) < threshold

2、三种梯度下降方法

（1）批量梯度下降
容易得到最优解，但每次都要考虑所有样本，速度很慢
（2）随机梯度下降
每次找一个样本，迭代速度快，但不一定每次都朝着收敛方向
（3）小批量梯度下降
每次更新选择一小部分数据来算

2.1 洗牌

import numpy.random
# import numpy.random# 对数据进行洗牌，shuffle函数是数据传进去就可以洗牌
def shuffleData(data):
    np.random.shuffle(data)
    cols = data.shape[1]
    X = data[:, 0:cols-1]
    y = data[:, cols-1:]
    return X, y

2.2 梯度下降求解

import time
# bathSize为1，则是随机梯度下降；为总的样本数，则是批量梯度下降；指定正1到总体之间，那就是小批量梯度下降
# thresh是阀值；alipha是学习率
def descent(data, theta, batchSize, stopType, thresh, alpha):
    #梯度下降求解
    
    init_time = time.time()
    i = 0 # 迭代次数
    k = 0 # batch
    X, y = shuffleData(data)
    grad = np.zeros(theta.shape) # 计算的梯度
    costs = [cost(X, y, theta)] # 损失值

    
    while True:
        grad = gradient(X[k:k+batchSize], y[k:k+batchSize], theta)
        k += batchSize #取batch数量个数据
        if k >= n: 
            k = 0 
            X, y = shuffleData(data) #重新洗牌
        theta = theta - alpha*grad # 参数更新
        costs.append(cost(X, y, theta)) # 计算新的损失
        i += 1 

        if stopType == STOP_ITER:       value = i
        elif stopType == STOP_COST:     value = costs
        elif stopType == STOP_GRAD:     value = grad
        if stopCriterion(stopType, value, thresh): break
    
    return theta, i-1, costs, grad, time.time() - init_time

2.3 画图

根据不同参数画图,根据选择传进来的参数,如果数据传进来一个,则是随机梯度下降方式,如果传进来的是总体则是批量梯度下降;如果传进来的是小批量的数据则是小批量梯度下降。

def runExpe(data, theta, batchSize, stopType, thresh, alpha):
	#先对一个值进行初始化,然后进行求解
    theta, iter, costs, grad, dur = descent(data, theta, batchSize, stopType, thresh,alpha)
    name = "Original" if (data[:,1]>2).sum() > 1 else "Scaled"
    name += " data - learning rate: {} - ".format(alpha)
    if batchSize==n: strDescType = "Gradient" # 批量梯度
    elif batchSize==1:  strDescType = "Stochastic" # 随机梯度
    else: strDescType = "Mini-batch ({})".format(batchSize) # 小批量梯度
    name += strDescType + " descent - Stop: "
    if stopType == STOP_ITER: strStop = "{} iterations".format(thresh)
    elif stopType == STOP_COST: strStop = "costs change < {}".format(thresh)
    else: strStop = "gradient norm < {}".format(thresh)
    name += strStop
    print ("***{}\nTheta: {} - Iter: {} - Last cost: {:03.2f} - Duration: {:03.2f}s".format(
        name, theta, iter, costs[-1], dur))
    fig, ax = plt.subplots(figsize=(12,4))
    ax.plot(np.arange(len(costs)), costs, 'r')
    ax.set_xlabel('Iterations')
    ax.set_ylabel('Cost')
    ax.set_title(name.upper() + ' - Error vs. Iteration')
    return theta

2.4 对比不同的停止策略

当n值指定为100的时候,相当于整体对于梯度下降,因为我的数据样本就100个.
（1）设置迭代次数停止策略
传进来的数据是按照迭代次数进行停止的,指定迭代次数的参数是thresh=5000.学习率是alpha=0.000001

n=100
runExpe(orig_data, theta, n, STOP_ITER, thresh=5000, alpha=0.000001)

（2）设置损失值停止策略
设置阀值为0.000001，当梯度值小于阀值的时候，进行停止

runExpe(orig_data, theta, n, STOP_COST, thresh=0.000001, alpha=0.001)

(3)根据梯度变化停止
设定阈值 0.05

runExpe(orig_data, theta, n, STOP_GRAD, thresh=0.05, alpha=0.001)

2.5 对比不同的梯度下降方法

（1）随机梯度下降，只有一个样本

runExpe(orig_data, theta, 1, STOP_ITER, thresh=5000, alpha=0.001)

把学习率调小，速度快，但稳定性差，需要很小的学习率

runExpe(orig_data, theta, 1, STOP_ITER, thresh=15000, alpha=0.000002)

（2）小批量梯度下降

runExpe(orig_data, theta, 16, STOP_ITER, thresh=15000, alpha=0.001)

学习率比较小,但是也有一些问题,浮动也很大,这个问题应该怎么解决呢?不再把学习率调小一些,而是换一种方案。
对数据进行标准化：将数据按其属性(按列进行)减去其均值，然后除以其方差。最后得到的结果是，对每个属性/每列来说所有数据都聚集在0附近，方差值为1。

from sklearn import preprocessing as pp
scaled_data = orig_data.copy()
scaled_data[:, 1:3] = pp.scale(orig_data[:, 1:3])
n=100
runExpe(scaled_data, theta, n, STOP_ITER, thresh=5000, alpha=0.001)