吴恩达机器学习：LR回归与正则化

小白初学

(一）介绍Logistic回归算法模型

1.1 LR回归

LR回归，虽然这个算法从名字上来看，是回归算法，但其实际上是一个分类算法，学术界也叫它logit regression, maximum-entropy classification (MaxEnt)或者是the log-linear classifier。在机器学习算法中，有几十种分类器，LR回归是其中最常用的一个。

问题引入：

我们知道，线性回归的模型是求出输出特征向量Y和输入样本矩阵X之间的线性关系系数θ，满足Y=Xθ。此时我们的Y是连续的，所以是回归模型。如果我们想要Y是离散的话，怎么办呢？

一个可以想到的办法是，我们对于这个Y再做一次函数转换，变为g(Y)。如果我们令g(Y)的值在某个实数区间的时候是类别A，在另一个实数区间的时候是类别B，以此类推，就得到了一个分类模型。如果结果的类别只有两种，那么就是一个二元分类模型了。逻辑回归的出发点就是从这来的。

1. LR回归是在线性回归模型的基础上，使用$sigmoid$函数，将线性模型 $w^{T}x$的结果压缩到 [0,1] 之间，使其拥有概率意义。 其本质仍然是一个线性模型，实现相对简单。

2. LR回归属于概率性判别式模型，之所谓是概率性模型，是因为LR模型是有概率意义的；之所以是判别式模型，是因为LR回归并没有对数据的分布进行建模，也就是说，LR模型并不知道数据的具体分布，而是直接将判别函数，或者说是分类超平面求解了出来。

一般来说，分类算法都是求解$p(C_k|x)$，即对于一个新的样本，计算其条件概率 $p(C_k|x)$。这个可以看做是一个后验概率，其计算可以基于贝叶斯公式得到：$p(C_k|x)=\frac{p(x|C_k)p(C_k)}{{\sum }_{k=1}^{K}p(x|C_k)p(C_k)}$ ,其中$p(x|C_k)$是类条件概率密度，$p(C_k)$是类的概率先验。使用这种方法的模型，称为是生成模型，即：$p(C_k|x)$是由$p(x|C_k)$和$p(C_k)$生成的。

分类算法所得到的 $p(C_k|x)$ 可以将输入空间划分成许多不相交的区域，这些区域之间的分隔面被称为判别函数（也称为分类面），有了判别函数，就可以进行分类了，上面生成模型，最终也是为了得到判别函数。如果直接对判别函数进行求解，得到判别面，这种方法，就称为判别式法。LR就属于这种方法。

总结： 逻辑回归是解决二分类问题的，简单来说就是 线性回归+sigmoid函数。

1.2 解析数学原理推导

详细参考：Logistic Regression 的前世今生（理论篇）

1.2.1 Logistic Distribution

LR回归是在线性回归模型的基础上，使用$sigmoid$函数得到的,先介绍一下$sigmoid$函数:

首先，需要对 logistic distribution 进行说明，这个分布的概率密度函数和概率分布函数如下：记$\mu$是位置参数，而s是形状参数

logistic distribution 分布在不同的 $\mu$ 和 s 的情况下，其概率密度函数 p(x;$\mu$,s) 的图形:

logistic distribution 在不同的 $\mu$ 和 s 的情况下，其概率分布函数 p(x;$\mu$,s) 的图形：

由图见，
logistic distribution和高斯分布的密度函数长得差不多。特别注意 logistic distribution的概率密度函数自中心附近增长速度较快，而在两端的增长速度相对较慢。形状参数s的数值越小，则概率密度函数在中心附近增长越快。

当 $\mu =0,s=1$ 时， logistic distribution的概率分布函数（累计密度函数）就是我们常说的$sigmoid$函数：

且其导数为：

这是一个非常好的特性且在后面的推导中将会被用到。

1.2.1 Logistic Regression代价函数

前面说过，分类算法都是求解$p(C_k|x)$，而Logistic Regression模型，就是用当$\mu =0,s=1$时的 logistic distribution的概率分布函数：$sigmoid$函数，对$p(C_k|x)$进行建模，所得到的模型，对于二分类的逻辑斯蒂回归模型有：

易求得：

这里$p$又可以被称作$odds$比值比，又名机会比、优势比。

$p$ 就是我们在分类算法中的决策面，由于LR回归是一个线性分类算法，所以，此处采用线性模型：这里参数向量为 ω = { w 0 , w 1 , w 2 , . . . , w n } ， ϕ j ( x ) \omega=\{w_{0},w_{1},w_{2},...,w_{n}\}，\phi_{j}(x) ω={w0​,w1​,w2​,...,wn​}，ϕj​(x)是线性模型的基函数，基函数的数目为M个，如果定义${\varphi }_0(x)=1$的话:



这里的p即$h_{\omega }(x)$（Hypothesis）。

那么， Logistic Regression模型$p(C_k|x)$就可以重写为下面这个形式：

回顾下线性回归的代价函数，由于线性回归是连续的，所以可以使用模型误差的的平方和来定义代价函数。但是逻辑回归不是连续的，线性回归代价函数定义的经验就用不上了。不过我们可以用最大似然法来推导出代价函数。

对于一个二分类的数据集 { x n , t n } \{x_n,t_n\} {xn​,tn​}，这里 t n ∈ { 0 , 1 } ， ϕ n = ϕ ( x n ) , t = { t 1 , t 2 , . . . , } , y n = p ( C 1 ∣ x n ) t_n\in\{0,1\}，\phi_n=\phi(x_n),t=\{t_1,t_2,...,\},y_n=p(C_1|x_n) tn​∈{0,1}，ϕn​=ϕ(xn​),t={t1​,t2​,...,},yn​=p(C1​∣xn​)其最大似然估计为：

对数似然估计为：（通常我们会取其平均1/N)

关于w的梯度为：


得到梯度之后，那么就可以和线性回归模型一样，使用梯度下降法求解参数$w$:

梯度下降法实现相对简单，但是其收敛速度往往不尽人意，可以考虑使用随机梯度下降法来解决收敛速度的问题。但上面两种在最小值附近，都存在以一种曲折的慢速逼近方式来逼近最小点的问题。所以在LR回归的实际算法中，用到的是牛顿法，拟牛顿法（DFP、BFGS、L-BFGS）。

由于求解最优解的问题，其实是一个凸优化的问题，这些数值优化方法的区别仅仅在于选择什么方向走向最优解，而这个方向通常是优化函数在当前点的一阶导数（梯度）或者二阶导数（海森Hessian矩阵）决定的。比如梯度下降法用的就是一阶导数，而牛顿法和拟牛顿法用的就是二阶导数。

1.3 逻辑回归对比线性回归

虽然逻辑回归能够用于分类，不过其本质还是线性回归。它仅在线性回归的基础上，在特征到结果的映射中加入了一层$sigmoid$函数（非线性）映射，即先把特征线性求和，然后使用$sigmoid$函数来预测。然而，正是这个简单的逻辑函数，使得逻辑回归模型成为了机器学习领域一颗耀眼的明星。

这主要是由于线性回归在整个实数域内敏感度一致，而分类范围，只需要在[0,1]之内。而逻辑回归就是一种减小预测范围，将预测值限定为[0,1]间的一种回归模型。逻辑曲线在z=0时，十分敏感，在$z>>0$或$z<<0$处，都不敏感，将预测值限定为[0,1]。

线性回归和LR的梯度更新公式是一样的,如下：

从梯度更新视角来看,唯一不同的是假设函数 h θ ( x ( i ) ) ， y ∈ { 0 , 1 } h_\theta(x^{(i)})，y\in \{0,1\} hθ​(x(i))，y∈{0,1}，对于LR而言，$h_{\theta }(x^{(i)})\in [0,1]$，那么梯度更新的幅度就不会太大。而线性回归$h_{\theta }(x^{(i)})$在整个实数域上，即使已经分类正确的点，在梯度更新的过程中也会大大影响到其它数据的分类，就比如有一个正样本，其输出为10000，此时梯度更新的时候，这个样本就会大大影响负样本的分类。而对于LR而言，这种非常肯定的正样本并不会影响负样本的分类情况。

(二）正则化

逻辑回归也会面临过拟合问题，所以我们也要考虑正则化。常见的有L1正则化和L2正则化。逻辑回归的L1正则化的损失函数表达式如下所述，相比普通的逻辑回归代价函数，增加了L1的范数做作为惩罚，超参数α作为惩罚系数，调节惩罚项的大小。

2.1 欠拟合 & 过度拟合

线性回归算法中：

逻辑回归算法中：

欠拟合是指忽略了训练集中某些数据，未能学习训练数据中的关系。
过拟合是指过分依赖训练数据，预测未知的数据的准确性较差。

解决方法：
（1）减少特征数量（减少特征会失去一些信息，即使特征选的很好）
• 可用人工选择要保留的特征。
• 模型选择算法。

（2）正则化（特征较多时比较有效）
• 保留所有特征，但减少 θ 的大小。

2.2 正则化

方法：

正则化是结构风险最小化策略的实现，是在经验风险上加一个正则化项或惩罚项。正则化项一般是模型复杂度的单调递增函数，模型越复杂，正则化项就越大。

线性回归模型的正则化：
代价函数变为：
$$J(\theta )=\frac{1}{2m}\sum _{i=1}^m(h\theta (x^{(i)})-y^{(i)}{)}^2+\lambda \sum _{j=1}^n{\theta }_j^2$$

公式前半部分是线性回归模型的代价函数，后半部分是正则化项。加上这个部分有两个目的，第一是维持对训练样本的拟合，第二是避免对训练样本的过度拟合。从数学方面来看，代价函数增加了一个正则化项后，代价函数不再仅仅由预测值和真实值的误差所决定，还和参数 $\theta$的大小有关。有了这个限制之后，要实现代价函数最小的目的，参数 $\theta$就不能随便取值，通过调节参数 $\lambda$，控制正则化项，避免过度拟合的情况。

正则化后的梯度下降变换为：
$${\theta }_0={\theta }_0-\alpha \frac{1}{m}\sum _{i=1}^m(h\theta (x^{(i)})-y^{(i)})\ast x_0^{(i)}$$

$${\theta }_j={\theta }_j-\alpha \frac{1}{m}\sum _{i=1}^m(h\theta (x^{(i)})-y^{(i)})\ast x_0^{(i)}+\frac{\lambda }{m}{\theta }_j$$

其中，$j=1,2,...,n$

合并后的式子为:

$${\theta }_j={\theta }_j(1-\frac{\alpha \lambda }{m})-\alpha \frac{1}{m}\sum _{i=1}^m(h\theta (x^{(i)})-y^{(i)})\ast x_0^{(i)}+\frac{\lambda }{m}{\theta }_j$$

逻辑回归模型的正则化：

对逻辑回归模型的代价函数进行正则化，其方法也是在原来的代价函数的基础上加上正则化项：

$$J(\theta )=-\frac{1}{m}[\sum _{i=1}^m{y}^{(i)}\cdot log(h\theta (x^{(i)}))+(1-y^{(i)})\cdot log(1-h\theta (x^{(i)}))]+\frac{1}{2m}\sum _{j=1}^n{\theta }_j^2$$

正则化后的梯度下降公式为：

$${\theta }_j={\theta }_j(1-\frac{\alpha \lambda }{m})-\alpha \frac{1}{m}\sum _{i=1}^m(h\theta (x^{(i)})-y^{(i)})\ast x_0^{(i)}+\frac{\lambda }{m}{\theta }_j$$

其中，$j=1,2,...,n$ 。由于${\theta }_0$ ​没有参与正则化。逻辑回归和线性回归看起来形式是一样的，但其实它们的算法并不一样，因为两个式子的预测函数 $h_{\theta }(x)$ 是不一样的。

此外 $\lambda$是正则项系数。当它的值很大时，说明对模型的复杂度惩罚大，对拟合数据的惩罚小，这样就不会出现过拟合情况。 在训练数据中的偏差比较大，在未知数据上的方差比较小，可能会出现欠拟合的情况。
同样，如果它的值小时，说明比较注重对训练数据的拟合，使得在训练数据上的偏差比较小，在未知数据上的方差较大，从而导致过拟合。

(三）LR回归算法作业

要求，建立一个logistic回归模型来预测一名学生是否能合格进入大学，根据两次考试的结果来决定每个申请人的录取机会。有以前的申请人的历史数据， 可以用它作为逻辑回归的训练集。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import matplotlib.ticker as ticker
import test
plt.style.use('fivethirtyeight') #样式美化
alpha = 0.004
iterations = 150000

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

x_plot = np.linspace(-10., 10., 5000)
y_plot = sigmoid(x_plot)
plt.plot(x_plot, y_plot)
plt.title('Sigmoid')
plt.show()

代价函数和梯度的求解：

def Cost_function(theta, X, y):
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    first = np.multiply(-y, np.log(sigmoid(X * theta.T)))
    second = np.multiply((1 - y), np.log(1 - sigmoid(X * theta.T)))
    return np.sum(first - second) / (len(X))

def predict(theta, X):
    probability = sigmoid(X * theta.T)
    return [1 if x >= 0.5 else 0 for x in probability]

def draw():
    # 定义x y数据 x1 y1:未通过 x2 y2:通过
    x1 = []
    y1 = []
    x2 = []
    y2 = []

导入训练数据：

 train_data = open("D:\机器学习\Andrew Ng Machine Learning\mlclass-ex2\ex2data1.txt")
    lines = train_data.readlines()
    for line in lines:
        scores = line.split(",")
        # 去除标记后面的换行符
        isQualified = scores[2].replace("\n", "")
        # 根据标记将两次成绩放到对应的数组
        if isQualified == "0":
            x1.append(float(scores[0]))
            y1.append(float(scores[1]))
        else:
            x2.append(float(scores[0]))
            y2.append(float(scores[1]))

    # 设置样式参数,默认主题 darkgrid（灰色背景+白网格）,调色板 2色
    sns.set(context="notebook", style="darkgrid", palette=sns.color_palette("RdBu", 2), color_codes=False)

创建两个分数的散点图：

positive = data[data.admitted.isin(['1'])]  # 1
negetive = data[data.admitted.isin(['0'])]  # 0

fig, ax = plt.subplots(figsize=(6,5))
ax.scatter(positive['exam1'], positive['exam2'], c='b', label='Admitted')
ax.scatter(negetive['exam1'], negetive['exam2'], s=50, c='r', marker='x', label='Not Admitted')
# 设置图例显示在图的上方
box = ax.get_position()
ax.set_position([box.x0, box.y0, box.width , box.height* 0.8])
ax.legend(loc='center left', bbox_to_anchor=(0.2, 1.12),ncol=3)
# 设置横纵坐标名
ax.set_xlabel('Exam 1 Score')
ax.set_ylabel('Exam 2 Score')
plt.show()

显示：横纵坐标是两门考试分数。蓝色代表通过。

path = 'D:\机器学习\Andrew Ng Machine Learning\mlclass-ex2\ex2data1.txt'
#读取CSV（逗号分割）文件到 DataFrame
data = pd.read_csv(path, header=None, names=['Exam 1', 'Exam 2', 'Admitted'])

#查看数据集的前 5 行
data.head()
print(data.head())

#查看数据的分布情况
data.describe()
print(data.describe())

数据预处理：

# add ones column
data.insert(0, 'Ones', 1)
# set X (training data) and y (target variable)
cols = data.shape[1]
X = data.iloc[:,0:cols-1]
y = data.iloc[:,cols-1:cols]
# convert to matrices and initialize theta
X = np.array(X.values)
y = np.array(y.values)
theta = np.matrix(np.array([0,0,0]))
#X.shape, theta.shape, y.shape
print('X.shape:',X.shape)
print('theta.shape:',theta.shape)
print('y.shape:',y.shape)
#计算一下初始损失值大小
Cost_function(theta, X, y)
print('Init_Cost_function:',Cost_function(theta, X, y))

梯度下降函数：

def gradient(theta, X, y):
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    parameters = int(theta.ravel().shape[1])
    grad = np.zeros(parameters)
    #误差
    error = sigmoid(X * theta.T) - y
    for i in range(parameters):
        #点乘 元素相乘
        term = np.multiply(error, X[:,i])
        grad[i] = np.sum(term) / len(X)
    return grad

gradient(theta,X,y)
print('gradient:',gradient(theta,X,y))

#func为优化目标函数，x0为依赖的参数，fprime为梯度下降函数
import scipy.optimize as opt   #行速度通常远远超过梯度下降
result = opt.fmin_tnc(func=Cost_function, x0=theta, fprime=gradient, args=(X, y))
print('result:',result)

#result = opt.minimize(fun=cost, x0=theta, args=(X, y), method='Newton-CG', jac=gradient)
#print(result)

# 设置训练得到的模型对应的直线，即h(x)对应的直线
# 设置x的取值范围：[30, 110]步长为10
#X = np.arange(30, 110, 10)
#y = (-result.X[0] - result.X[1] * Xx) / result.X[2]
#plt.figure()
#plt.plot(X, y)

代入参数值，计算损失函数的值：

Cost_function(result[0], X, y)
print('Cost_function:',Cost_function(result[0], X, y))

预测准确率评估：

theta_min = np.matrix(result[0])
predictions = predict(theta_min, X)
print('predictions:',predictions)

correct = [1 if ((a == 1 and b == 1) or (a == 0 and b == 0)) else 0 for (a, b) in zip(predictions, y)]
#map()函数将correct转化为全部为int的列表
accuracy = sum(correct) / len(X)
print ('accuracy = {0}%'.format(accuracy))

# 绘图
draw()

实验结果：

sigmoid 函数图：
显示数据中前5行以及数据的分布情况：

绘制的散点图：

X，y，theta的维数：

计算得到初始的代价函数的值：

批量梯度下降：

梯度下降优化后的代价函数，寻找的参数 $\theta$值：

带入参数值计算得到的代价函数的值：

用训练集进行预测：

预测结果的准确率：

0.89 #预测精度达到了89%

寻找到的决策边界：

(四）LR回归正则化算法作业

设想你是工厂的生产主管，你有一些芯片在两次测试中的测试结果。对于这两次测试，你想决定是否芯片要被接受或抛弃。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

data2 = pd.read_csv('D:\机器学习\Andrew Ng Machine Learning\mlclass-ex2\ex2data2.txt', names=['Test 1', 'Test 2', 'Accepted'])
data2

画出散点图，观察是否呈线性：

def plot_data():
    positive = data2[data2['Accepted'].isin([1])] #将Accepted列中是1的数据提出来，单独呈一个表格
    negative = data2[data2['Accepted'].isin([0])] #同上，提出Accepted是0的表格
    #接下来画图
    fig, ax = plt.subplots(figsize=(8,5))
    ax.scatter(positive['Test 1'], positive['Test 2'], s=50, c='b', marker='o', label='Accepted')
    ax.scatter(negative['Test 1'], negative['Test 2'], s=50, c='r', marker='x', label='Rejected')
    ax.legend() #画出图例
    ax.set_xlabel('Test 1 Score')
    ax.set_ylabel('Test 2 Score')


plot_data()
plt.show() 

分析：观察上图发现，正负两类数据并没有线性的决策界限，因此直接用logistic回归算法在这个数据集上并不能表现良好。因此需要进行正则化。

定义一个特征映射函数：

def feature_mapping(x1, x2, power):
    data = {}    #定义一个字典，从而在之后可以直接添加元素的名称和内容
    for i in np.arange(power + 1):    #i=n时，就是n次幂的所有组合。一共最多power次幂
        for p in np.arange(i + 1):
            data["f{}{}".format(i - p, p)] = np.power(x1, i - p) * np.power(x2, p)  #format是格式化函数。
    return pd.DataFrame(data)   

初始化数据：

x1 = data2['Test 1'].values  #定义x1，x2的数据
x2 = data2['Test 2'].values
_data2 = feature_mapping(x1, x2, power=6)  #得到特征映射后的二维表
_data2.head()
X = _data2.values   #将二维表转化为array格式
y = data2['Accepted'].values #初始化y
theta = np.zeros(X.shape[1])
X.shape, y.shape, theta.shape 

定义代价函数：正则化的代价函数，解决过拟合现象

def cost(theta, X, y):    #这是未正则化的代价函数 
    first = (-y) * np.log(sigmoid(X @ theta)) #注意array中有三种乘法
    second = (1 - y)*np.log(1 - sigmoid(X @ theta))
    return np.mean(first - second)

def sigmoid(z):  # logist函数
    return 1 / (1 + np.exp(- z))

def costReg(theta, X, y, l=1):   #这是正则化之后的代价函数
    _theta = theta[1: ]  #不惩罚第一项
    reg = (l / (2 * len(X))) *(_theta @ _theta)  
    return cost(theta, X, y) + reg  

定义批量的梯度函数：同样是正则化后的那个导数

def gradient(theta, X, y):   #未正则化 
    return (X.T @ (sigmoid(X @ theta) - y))/len(X)
    
def gradientReg(theta, X, y, l=1):  #正则化后 
    reg = (1 / len(X)) * theta
    reg[0] = 0    #同样不惩罚第一个θ
    return gradient(theta, X, y) + reg
    gradientReg(theta, X, y, 1)   
#通过高级优化来进行梯度下降（计算速度很快，且不需要人为设定α）
import scipy.optimize as opt
result2 = opt.fmin_tnc(func=costReg, x0=theta, fprime=gradientReg, args=(X, y, 2))   

#判断拟合的模型结果的精确性
final_theta = result2[0]
def predict(theta, X):
    probability = sigmoid(X@theta)
    return [1 if x >= 0.5 else 0 for x in probability] 
predictions = predict(final_theta, X)
correct = [1 if a==b else 0 for (a, b) in zip(predictions, y)]  #通过模型预测的结果与真实结果相比
accuracy = sum(correct) / len(correct)

画出决策边界：

1. 初始化网格点
2. 所有网格点变成高次幂的形式、再×theta。这个式子命名为z，z等于0时，y=0.5。contour函数是默认z=0

x = np.linspace(-1, 1.5, 250)  #x维度是250
xx, yy = np.meshgrid(x, x)  #xx维度是250,250  xx.ravel()的维度是62500，ravel即多维数组转化为一维数组,后面featuremapping的参数是要一维的。
z = feature_mapping(xx.ravel(), yy.ravel(), 6).values  #z的维度是62500,28
z = z @ final_theta  #z的维度是62500
z = z.reshape(xx.shape)  #z的维度是250,250
plot_data()  # 开始定义的函数，执行后画出散点图.
plt.contour(xx, yy, z, 0)  #这里z是一个关于xx、yy的高次幂的函数，默认z=0画图
plt.ylim(-.8, 1.2)
plt.show()   

预测结果的准确率：

0.8050847457627118

寻找到的决策边界：
(用本题的最后画图的contour(xx, yy, z, 0)方法来画第一题的图，结果和第一题完全一样。)