练习二 - 逻辑回归(Logistic regression)

Detection
293 阅读12分钟

Logistic regression

一、引言

本次练习中，我们将会实现逻辑回归并将其运用于两个不同的数据集。在开始本次联系之前，需要导入这两个训练集(ex2data1.txt和ex2data2.txt)，都可从GitHub上获取。

二、逻辑回归

本节中，我们将会实现一个可以用来预测一个学生是否会被大学拟录取的逻辑回归模型。

1. 问题背景

假如你是一个大学部门的管理人员并且你想要基于每个申请学生的两次测验成绩来决定其是否被录取。当然，你拥有被录取学生的历史数据，这些数据可以作为我们逻辑回归的训练集。每个训练样本包括了申请者两次测验的分数以及录取结果。

我们所要完成的任务就是建立一个分类模型，基于两次测验成绩来估算每个申请者被录取的可能性。

2. 具体实现

A. 导入数据

让我们导入数据并检验

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

path = 'ex2data1.txt'
data = pd.read_csv(path, header=None, names=['Exam 1', 'Exam 2', 'Admitted'])
# 查看前五个样本数据
data.head()
Exam 1Exam 2Admitted
034.62366078.0246930
130.28671143.8949980
235.84740972.9021980
360.18259986.3085521
479.03273675.3443761

让我们创建两测验分数的散点图，并使用颜色编码来可视化，如果样本是正的（被录取）或负的（未被录取）。

# 利用切片分离正负样本
positive = data[data['Admitted'].isin([1])]
negative = data[data['Admitted'].isin([0])]
​
fig, ax = plt.subplots(figsize=(12,8))
# 绘出正样本点的散点图
ax.scatter(positive['Exam 1'], positive['Exam 2'], s=50, c='b', marker='o', label='Admitted')
# 绘出负样本点的散点图
ax.scatter(negative['Exam 1'], negative['Exam 2'], s=50, c='r', marker='x', label='Not Admitted')
ax.legend()
ax.set_xlabel('Exam 1 Score')
ax.set_ylabel('Exam 2 Score')
plt.show()

download.png

从上图我们可以看出两种分类有着清晰的决策边界。现在我们需要实现逻辑回归，那样就可以训练一个模型来预测结果。方程实现在下面的代码示例在"exercises" 文件夹的 "ex2.pdf" 中。

B. Sigmoid函数

g 代表一个常用的逻辑函数（logistic function）为S形函数（Sigmoid function），公式为：

g(z)=11+e−zg(z) = \frac{1}{1+e^{-z}}

合起来，我们得到逻辑回归模型的假设函数：

hθ(x)=g(θTx)=11+e−θTXh_\theta(x) = g(\theta^Tx) =\frac{1}{1 +e^{ -\theta^T X}}

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

试一试这个函数，看是否能正确得出结果

nums = np.arange(-10, 10, step=1)
​
fig, ax = plt.subplots(figsize=(12,8))
# 用20段线段构成
ax.plot(nums, sigmoid(nums), 'r')
plt.show()

download.png

效果看起来很不错。接下来我们需要定义代价函数对预测结果进行评估。代价函数如下：

C. 代价函数

J(θ)=1m∑i=1m[−y(i)log⁡(hθ(x(i)))−(1−y(i))log⁡(1−hθ(x(i)))]J(\theta)=\frac{1}{m}\sum_{i=1}^{m}\left[-y^{(i)} \log \left(h_{\theta}\left(x^{(i)}\right)\right)-\left(1-y^{(i)}\right) \log \left(1-h_{\theta}\left(x^{(i)}\right)\right)\right]

同时我们也可以得到代价函数对每个参数的梯度(偏导数)：

∂J(θ)∂θj=1m∑i=1m(hθ(x(i))−y(i))xj(i)\frac{\partial J(\theta)}{\partial \theta_{j}}=\frac{1}{m} \sum_{i=1}^{m}\left(h_{\theta}\left(x^{(i)}\right)-y^{(i)}\right) x_{j}^{(i)}

def cost(theta, X, y):
    theta = np.matrix(theta)
    # 转换为Numpy矩阵便于矩阵运算
    X = np.matrix(X)
    y = np.matrix(y)
    # 注意此处的X和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))

现在，我们要做一些设置，和我们在练习1在线性回归的练习很相似。

# 在样本集中添加一列全为一的数据便于矩阵乘法
data.insert(0, 'Ones', 1)
​
# 分离特征向量和对应标签
cols = data.shape[1]
X = data.iloc[:,0:cols-1]
y = data.iloc[:,cols-1:cols]
​
# 转换X和y的数据类型
X = np.array(X.values)
y = np.array(y.values)
​
# 对参数矩阵进行初始化，原特征向量维数为2，添加一列后变成3，那么参数矩阵的列数就为 3
theta = np.zeros(3)

检查一下各参数维度：

theta

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

X.shape, theta.shape, y.shape

((100, 3), (3,), (100, 1))

让我们计算初始化参数的代价函数(theta为0)。

cost(theta, X, y)

0.6931471805599453

由于 Θ = 0,那么 z=θTx=0,g(z)=Sigmoid(0)=0.5z = \theta^Tx = 0,g(z) = Sigmoid(0) = 0.5 。所以我们计算的是：J(0)=1100∑i=1100[−y(i)log(0.5)−(1−y(i))log(0.5)]J(0) = \frac{1}{100} \sum^{100}_{i = 1 }[-y^{(i)}log(0.5) - (1-y^{(i)})log(0.5)]

# 化简可得
np.log(2)

0.6931471805599453

看起来不错，接下来，我们需要一个函数来计算我们的训练数据、标签和一些参数thata的梯度。

D. 梯度下降(gradient descent)

之前提到过代价函数对各个参数的偏导数： ∂J(θ)∂θj=1m∑i=1m(hθ(x(i))−y(i))xj(i)\frac{\partial J(\theta)}{\partial \theta_{j}}=\frac{1}{m} \sum_{i=1}^{m}\left(h_{\theta}\left(x^{(i)}\right)-y^{(i)}\right) x_{j}^{(i)}

这些导数看起来和线性回归的梯度十分相似，但是由于线性回归和逻辑回归定义的假设函数形式不同(线性回归用直线你和，逻辑回归使用Sigmoid函数预测概率)，所以公式内涵不同

  • 注意：
  1. 这是批量梯度下降
  2. 转化为向量进行计算：1mXT(Sigmoid(Xθ)−y)\frac{1}{m} X^T( Sigmoid(X\theta) - y )
def gradient(theta, X, y):
    # 先进行数据转换便于矩阵乘法
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    
    # 定义一个存储矩阵用于存储梯度，是一个和theta矩阵分量个数相同的一维矩阵
    parameters = int(theta.ravel().shape[1])
    grad = np.zeros(parameters)
    
    # 计算预测
    error = sigmoid(X * theta.T) - y
    
    for j in range(parameters):
        term = np.multiply(error, X[:,j])
        grad[j] = np.sum(term) / len(X)
    # 返回梯度行向量
    return grad

注意，我们实际上没有在这个函数中执行梯度下降，我们仅仅在计算一个梯度步长。在练习中，一个称为“fminunc”的Octave函数是用来优化函数来计算成本和梯度参数。由于我们使用Python，我们可以用SciPy的“optimize”命名空间来做同样的事情。

我们看看用我们的数据和初始参数为0的梯度下降法的结果。

gradient(theta, X, y)

array([ -0.1       , -12.00921659, -11.26284221])

现在可以用SciPy's truncated newton（TNC）实现寻找最优参数。

  • 最常使用的参数：
  1. func：优化的目标函数
  2. x0：初值
  3. fprime：提供优化函数func的梯度函数，不然优化函数func必须返回函数值和梯度，或者设置approx_grad=True
  4. approx_grad :如果设置为True，会给出近似梯度
  5. args：元组，是传递给优化函数的参数
  • 返回：
  1. x ： 数组，返回的优化问题目标值
  2. nfeval ： 整数，function evaluations的数目在进行优化的时候，每当目标优化函数被调用一次，就算一个function evaluation。在一次迭代过程中会有多次function evaluation。这个参数不等同于迭代次数，而往往大于迭代次数。
  3. rc ： int,Return code, see below
import scipy.optimize as opt
result = opt.fmin_tnc(func=cost, x0=theta, fprime=gradient, args=(X, y))
result

(array([-25.16131863,   0.20623159,   0.20147149]), 36, 0)

得到最优参数以后，我们再来计算一下代价函数值：

cost(result[0],X,y)

0.20349770158947458

E. 得出分类结果

接下来，我们需要编写一个函数，用我们所学的参数theta来为数据集X输出预测。然后，我们可以使用这个函数来给我们的分类器的训练精度打分。 逻辑回归模型的假设函数：

hθ(x)=11+e−θTXh_\theta(x) = \frac{1}{1+e^{-\theta^TX }}

当 ℎ_𝜃 大于等于0.5时，预测 y=1

当 ℎ_𝜃 小于0.5时，预测 y=0 。

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

# 得到使代价函数最小的Θ
theta_min = np.matrix(result[0])
# 直接最原训练集中的所有数据进行预测，得到结果向量和标签向量y尺寸一致
predictions = predict(theta_min, X)
# 利用列表生成式逐个得出预测准确与否
correct = [1 if ((a == 1 and b == 1) or (a == 0 and b == 0)) else 0 for (a, b) in zip(predictions, y)]
# 计算预测精度
accuracy = (sum(map(int, correct)) % len(correct))
print ('accuracy = {0}%'.format(accuracy))

accuracy = 89%

我们的逻辑回归分类器预测正确，如果一个学生被录取或没有录取，达到89%的精确度。不坏！记住，这是训练集的准确性。我们没有保持住了设置或使用交叉验证得到的真实逼近，所以这个数字有可能高于其真实值（这个话题将在以后说明）。

F. 如何绘制决策边界？

首先我们要思考决策边界从何而来？我们在逻辑回归中，以0.5作为界限，利用预测出来的正负标签的概率直接得出结果，即：

ypredict={1,hθ(z)≥0.50,hθ(z)<0.5y_{predict} = \begin{cases} 1 ,h_\theta(z) \ge 0.5\\ 0,h_\theta(z) < 0.5 \end{cases}

又z=θTX，hθ(z)=Sigmoid(z)z = \theta^TX， h_\theta(z) = Sigmoid(z)，所以0.5这个界限转化为特征向量关系就是 z=θTX=0z = \theta^TX = 0

带入此问题，也就是：θ0+θ1x1+θ2x2=0 \theta_0 + \theta_1 x_1 + \theta_2 x_2 = 0

由于Θ已知，所以我们的函数可以表示为：x2=−1θ2(θ0+θ1x1)x_2 = -\frac{1}{\theta_2} (\theta_0 + \theta_1 x_1)

# 注意参数的尺寸(1,3)
theta_min

matrix([[-25.16131863,   0.20623159,   0.20147149]])

# 利用切片分离正负样本
positive = data[data['Admitted'].isin([1])]
negative = data[data['Admitted'].isin([0])]

# 得出直线上两点
point1 = (0, -theta_min[0,0]/theta_min[0,2])
point2 = (-theta_min[0,0]/theta_min[0,1],0)

fig, ax = plt.subplots(figsize=(12,8))
# 绘出正样本点的散点图
ax.scatter(positive['Exam 1'], positive['Exam 2'], s=50, c='b', marker='o', label='Admitted')
# 绘出负样本点的散点图
ax.scatter(negative['Exam 1'], negative['Exam 2'], s=50, c='r', marker='x', label='Not Admitted')

# 绘出决策边界
ax.plot(point1,point2, c = 'k', label = 'Decision boundary')

ax.legend()
ax.set_xlabel('Exam 1 Score')
ax.set_ylabel('Exam 2 Score')
plt.show()

download.png

三、正则化逻辑回归

在训练的第二部分，我们将要通过加入正则项提升逻辑回归算法。如果你对正则化有点眼生，或者喜欢这一节的方程的背景，请参考在"exercises"文件夹中的"ex2.pdf"。简而言之，正则化是代价函数中的一个术语，它使算法更倾向于“更简单”的模型（其实也就是对高维x赋予更小的系数）。正则化助于减少过拟合，提高模型的泛化能力。这样，我们开始吧。

1. 问题背景

设想你是工厂的生产主管，你有一些芯片在两次测试中的测试结果。对于这两次测试，你想决定是否芯片要被接受或抛弃。为了帮助你做出艰难的决定，你拥有过去芯片的测试数据集，从其中你可以构建一个逻辑回归模型。

和第一部分类似，从数据可视化开始吧！

2. 导入数据

path =  'ex2data2.txt'
data2 = pd.read_csv(path, header=None, names=['Test 1', 'Test 2', 'Accepted'])
data2.head()
Test 1Test 2Accepted
00.0512670.699561
1-0.0927420.684941
2-0.2137100.692251
3-0.3750000.502191
4-0.5132500.465641
positive = data2[data2['Accepted'].isin([1])]
negative = data2[data2['Accepted'].isin([0])]

fig, ax = plt.subplots(figsize=(12,8))
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')
plt.show()

download.png

哇，这个数据看起来可比前一次的复杂得多。可以看出，其中没有可以用来分离两个分类的线性决策界限。一个方法是用像逻辑回归这样的线性技术来构造从原始特征的多项式中得到的特征。让我们通过创建一组多项式特征入手吧。

degree = 5
# 分离特征向量的分量
x1 = data2['Test 1']
x2 = data2['Test 2']

# 为原样本集插入一列全一，便于矩阵乘法
data2.insert(3, 'Ones', 1)

# 生成多项式函数f值
for i in range(1, degree):
    for j in range(0, i):
        data2['F' + str(i) + str(j)] = np.power(x1, i-j) * np.power(x2, j)

# 删去原来的
data2.drop('Test 1', axis=1, inplace=True)
data2.drop('Test 2', axis=1, inplace=True)

data2.head()
AcceptedOnesF10F20F21F30F31F32F40F41F42F43
0110.0512670.0026280.0358640.0001350.0018390.0250890.0000070.0000940.0012860.017551
111-0.0927420.008601-0.063523-0.0007980.005891-0.0435090.000074-0.0005460.004035-0.029801
211-0.2137100.045672-0.147941-0.0097610.031616-0.1024120.002086-0.0067570.021886-0.070895
311-0.3750000.140625-0.188321-0.0527340.070620-0.0945730.019775-0.0264830.035465-0.047494
411-0.5132500.263426-0.238990-0.1352030.122661-0.1112830.069393-0.0629560.057116-0.051818

注意我们对原始数据集特征结果的转换：我们将原特征转化为核函数：

F10=x11×x20F20=x12×x20F21=x11×x21F30=x13×x20F31=x12×x21F32=x11×x22F40=x14×x20F41=x13×x21F42=x12×x22F43=x11×x23F_{10} = x_1^1 \times x_2 ^0 \\ F_{20} = x_1^2 \times x_2 ^0 \\ F_{21} = x_1^1 \times x_2 ^1 \\ F_{30} = x_1^3 \times x_2 ^0 \\ F_{31} = x_1^2 \times x_2 ^1 \\ F_{32} = x_1^1 \times x_2 ^2 \\ F_{40} = x_1^4 \times x_2 ^0 \\ F_{41} = x_1^3 \times x_2 ^1 \\ F_{42} = x_1^2 \times x_2 ^2 \\ F_{43} = x_1^1 \times x_2 ^3 \\

现在，我们需要修改第1部分的代价和梯度函数，包括正则化项。首先是代价函数：

3. 正则化代价函数

J(θ)=1m∑i=1m[−y(i)log⁡(hθ(x(i)))−(1−y(i))log⁡(1−hθ(x(i)))]+λ2m∑j=1nθj2J\left( \theta \right)=\frac{1}{m}\sum\limits_{i=1}^{m}{[-{{y}^{(i)}}\log \left( {{h}_{\theta }}\left( {{x}^{(i)}} \right) \right)-\left( 1-{{y}^{(i)}} \right)\log \left( 1-{{h}_{\theta }}\left( {{x}^{(i)}} \right) \right)]}+\frac{\lambda }{2m}\sum\limits_{j=1}^{n}{\theta _{j}^{2}}

# 此处的正则化系数也可以看作是学习率
def costReg(theta, X, y, Lambda):
    # 先进行数据转换
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    
    # 矩阵与矩阵要使用矩阵乘法：np.multiply
    first = np.multiply(-y , np.log(sigmoid(X * theta.T)))
    second =np.multiply( -(1 - y) , np.log(1 - sigmoid(X * theta.T)))
    
    
    # 注意我们此处的θj是不包含 j = 0的，所以需要进行切片
    regula = Lambda/(2* len(X)) * np.sum(np.power(theta[:, 1:theta.shape[1]],2))
    return np.sum(first + second) / len(X) + regula

请注意等式中的"reggula" 项。还注意到另外的一个“学习率（lambda）”参数。这是一种超参数，用来控制正则化项。现在我们需要添加正则化梯度函数：

如果我们要使用梯度下降法令这个代价函数最小化，因为我们未对 𝜃_0 进行正则化，所以梯度下降算法将分两种情形：

Repeatuntilconvergence{θ0:=θ0−a1m∑i=1m[hθ(x(i))−y(i)]x0(i)θj:=θj−a1m∑i=1m[hθ(x(i))−y(i)]xj(i)+λmθj} Repeat \begin{array}{l} Repeat until convergence\{\\ \theta_{0}:=\theta_{0}-a \frac{1}{m} \sum_{i=1}^{m}\left[h_{\theta}\left(x^{(i)}\right)-y^{(i)}\right] x_{0}^{(i)} \\ \theta_{j}:=\theta_{j}-a \frac{1}{m} \sum_{i=1}^{m}\left[h_{\theta}\left(x^{(i)}\right)-y^{(i)}\right] x_{j}^{(i)}+\frac{\lambda}{m} \theta_{j} \\ \} \\ \text { Repeat } \end{array}

对上面的算法中 j=1,2,...,n 时的更新式子进行调整可得： θj:=θj(1−aλm)−a1m∑i=1m(hθ(x(i))−y(i))xj(i){{\theta }_{j}}:={{\theta }_{j}}(1-a\frac{\lambda }{m})-a\frac{1}{m}\sum\limits_{i=1}^{m}{({{h}_{\theta }}\left( {{x}^{(i)}} \right)-{{y}^{(i)}})x_{j}^{(i)}}

def gradientReg(theta, X, y, learningRate):
    # 修改数据类型
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    
    # 得到theta矩阵中的分量个数，定义一个空梯度向量，尺寸为(1, 2*len(theta))
    parameters = int(theta.ravel().shape[1])
    grad = np.zeros(parameters)
    
    # 求出预测误差向量
    error = sigmoid(X * theta.T) - y
    
    # 使用修改过后的梯度计算方式
    for j in range(parameters):
        # 行向量与X矩阵的第j列分量相乘
        term = np.multiply(error, X[:,j])
        
        if (j == 0):
            grad[j] = np.sum(term) / len(X)
        else:
            grad[j] = (np.sum(term) / len(X)) + ((learningRate / len(X)) * theta[:,j])
    
    return grad

就像在第一部分中做的一样，初始化变量。

# 设定特征矩阵X和标签向量y，注意此处我们之前将标签向量移动到了data2的第一列( i = 0)且
# 此处的特征向量并非原始特征向量而是利用核函数简化后的特征向量
cols = data2.shape[1]
X2 = data2.iloc[:,1:cols]
y2 = data2.iloc[:,0:1]

# 转变数据类型
X2 = np.array(X2.values)
y2 = np.array(y2.values)
# theta是一个行向量
theta2 = np.zeros(11)

theta2

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

查看一下数据尺寸是否符合题意

X2.shape,y2.shape,theta2.shape

((118, 11), (118, 1), (11,))

让我们初始学习率(正则化系数)到一个合理值。如果后面有必要的话（即如果惩罚太强或不够强）,我们可以之后再捣腾此值(抑制过拟合倾向)。

learningRate = 1

现在，让我们尝试调用新正则化函数，依旧将的默认为0的theta传入，以确保计算工作正常。

# 初始正则化代价函数值，lambda = 0 = theta,所以结果仍为 log2
costReg(theta2, X2, y2, learningRate)

0.6931471805599454

# 首次计算梯度值
gradientReg(theta2, X2, y2, learningRate)

array([0.00847458, 0.01878809, 0.05034464, 0.01150133, 0.01835599,
       0.00732393, 0.00819244, 0.03934862, 0.00223924, 0.01286005,
       0.00309594])

import scipy.optimize as opt

result2 = opt.fmin_tnc(func=costReg, x0=theta2, fprime=gradientReg, args=(X2, y2, learningRate))
result2

(array([ 0.53010247,  0.29075567, -1.60725764, -0.58213819,  0.01781027,
        -0.21329507, -0.40024142, -1.3714414 ,  0.02264304, -0.95033581,
         0.0344085 ]),
 22,
 1)

theta_min = np.matrix(result2[0])

# 将逻辑回归模型返回的概率转换为预测结果
predictions = predict(theta_min, X2)
correct = [1 if ((a == 1 and b == 1) or (a == 0 and b == 0)) else 0 for (a, b) in zip(predictions, y2)]

accuracy = (sum(map(int, correct)) / len(correct))
print ('accuracy = {0}%'.format(accuracy))

accuracy = 0.6610169491525424%

虽然我们实现了这些算法，值得注意的是，我们还可以使用高级Python库,如scikit-learn来解决这个问题。

# 调用sklearn的线性回归包
from sklearn import linear_model
# 乘法l2，学习率(正则化系数)初值为 1.0
model = linear_model.LogisticRegression(penalty='l2', C=1.0)
# 拟合
model.fit(X2, y2.ravel())

LogisticRegression()

model.score(X2, y2)

0.6610169491525424
avatar
文章
阅读
粉丝