逻辑回归练习
逻辑回归
练习1 无正则化逻辑回归
在训练的初始阶段，我们将要构建一个逻辑回归模型来预测，某个学生是否被大学录取。设想你是大学相关部分的管理者，想通过申请学生两次测试的评分，来决定他们是否被录取。现在你拥有之前申请学生的可以用于训练逻辑回归的训练样本集。对于每一个训练样本，你有他们两次测试的评分和最后是被录取的结果。为了完成这个预测任务，我们准备构建一个可以基于两次测试评分来评估录取可能性的分类模型。
可视化数据
导入需要用到的库
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
# 读取数据
path = "D:\gitcode\Coursera-ML-AndrewNg-Notes-master\Coursera-ML-AndrewNg-Notes-master\code\ex2-logistic regression\ex2data1.txt"
data = pd.read_csv(path, names=["exam1", "exam2", "admitted"])
# print(data.head())
让我们创建两个分数的散点图，并使用颜色编码来可视化，如果样本为1的（被接纳）或0的（未被接纳）。
# 可视化数据
positive = data[data['admitted'].isin([1])]
negative = data[data['admitted'].isin([0])]
fig, ax = plt.subplots(figsize=(6, 5))
ax.scatter(positive['exam1'], positive['exam2'], c='b', label='admitted')
ax.scatter(negative['exam1'], negative['exam2'], c='r', label='not admitted')
plt.xlabel("exam1 score")
plt.ylabel("exam2 score")
plt.legend()
plt.show()
sigmoid 函数
g 代表一个常用的逻辑函数（logistic function）为S形函数（Sigmoid function），公式为： \[g\left( z \right)=\frac{1}{1+{{e}^{-z}}}\]
得到逻辑回归模型的假设函数：
\[{{h}_{\theta }}\left( x \right)=\frac{1}{1+{{e}^{-X{{\theta }^{T}}}}}\]
定义sigmoid函数
# sigmoid函数
def sigmoid(z):
return 1 / (1 + np.exp(-z))
代价(损失)函数
现在，我们需要编写代价函数来评估结果。
代价函数：
\(J\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)]}\)
在用代码实现时，需要注意，theta维1维向量，使用numpy进行矩阵相乘时，需要是numpy数组，否则可能会存在问题
# 定义损失函数
def computeCost(theta, X, y):
'''
:param theta: 1维向量 n个特征
:param X: numpy数组 m行n列
:param y: numpy数组 m行1列
:return: 损失数值
'''
theta = np.reshape(theta, (1, theta.shape[0]), order='C')
m = X.shape[0]
first = (-y) * np.log(sigmoid(X @ theta.T))
second = (1 - y) * np.log(1 - sigmoid(X @ theta.T))
return np.sum(first - second) / m
梯度函数
- 这是梯度公式（batch gradient descent）：
# 计算梯度
def gradient(theta, X, y):
'''
:param theta: n维向量 n个特征
:param X: numpy数组 m行n列
:param y: numpy数组 m行n列
:return: n维梯度 n个特征(n,)
'''
theta = np.reshape(theta, (1, theta.shape[0]), order='C')
m = X.shape[0]
return (((sigmoid(X @ theta.T) - y).T @ X) / m).flatten()
在实现时，使用矩阵乘法进行得到梯度矩阵，最终返回得到一个n维向量，因为在后续使用tnc算法进行优化时，要求optimize的函数返回值为向量。此外，注意，我们实际上没有在这个函数中执行梯度下降，我们仅仅在计算一个梯度步长。在练习中，一个称为“fminunc”的Octave函数是用来优化函数来计算成本和梯度参数。由于我们使用Python，我们可以用SciPy的“optimize”命名空间来做同样的事情。
初始化变量
在进行TNC寻找最优参数之前，我们还需要进行变量的初始化，首先我们需要为数据都添加一列全为1的值，这代表偏置项的特征值，其次，我们还需要将dataframe转化成numpy数组，才能进行矩阵运算，否则，可能会出现不该出现的问题。
# 初始化变量
data.insert(0, 'Ones', 1)
X = data.drop(columns=['admitted']).values
y = data['admitted'].values.reshape(-1, 1)
theta = np.zeros(X.shape[1]).flatten()
SciPy's truncated newton（TNC）实现寻找最优参数
我们可以利用scipy.optimize提供的fmin_tnc进行寻找最优参数
# 利用tcn进行最优化拟合
# 方法一直接利用opt.fmin_tnc进行优化
import scipy.optimize as opt
result = opt.fmin_tnc(func=computeCost, x0=theta, fprime=gradient, args=(X, y))
theta = np.reshape(result[0], (1, 3), order='C')
也可以利用opt.minimize进行优化，设置优化方法为TNC
# 方法二利用opt.minimize进行优化,设置方法为TNC
result = opt.minimize(fun=computeCost, x0=theta, jac=gradient, args=(X, y), method='TNC')
theta = np.reshape(result['x'], (1,3), order='C')
print(result)
但是，注意，两种方法得到的结果的类型不同，可以将结果打印出来后取出所需要的最优参数
预测函数
接下来，我们需要编写一个函数，用我们所学的参数theta来为数据集X输出预测。然后，我们可以使用这个函数来给我们的分类器的训练精度打分。
逻辑回归模型的假设函数：
\[{{h}_{\theta }}\left( x \right)=\frac{1}{1+{{e}^{-X{{\theta }^{T}}}}}\]
当\({{h}_{\theta }}\)大于等于0.5时，预测y=1，
当\({{h}_{\theta }}\)小于0.5时，预测 y=0。
# 预测函数
def predict(theta, X):
'''
:param theta: 所学的参数theta n维向量
:param X: 所需拟合的参数
:return: 类别列表
'''
probability = sigmoid(X @ theta.T)
# print(probability)
return [1.0 if x >= 0.5 else 0.0 for x in probability]
将我们所需要预测的数据值做为X，theta所我们学习出来的参数，根据逻辑回归模型的假设函数进行返回类别列表。
我们使用训练集的输入X进行预测返回预测的类别列表，之后判断其是否和预期值一样，进而计算出训练精确度。
# 计算预测精度
prediction = predict(theta, X)
num = 0
for i in range(0, len(prediction)):
if (prediction[i] == y[i][0]):
num = num + 1
acc = num / len(X)
print(acc)
在这里我们得到无论是使用opt.minimize选择优化方法为tnc或者是直接使用opt.fmin_tnc得到的逻辑回归分类器预测正确，如果一个学生被录取或没有录取，达到89%的精确度。
绘制决策界限
为了更好展示逻辑回归所得到的决策边界的效果，我们进行绘制界限图，利用等高线图的特点，取在exam1范围内的网格点，利用矩阵相乘进行计算(\(X{{\theta }^{T}}\))后得到许多值，但是我们只需要等高线等于0的线，因为其对应的为我们的分界线。
当然，我们也可以直接通过exam1的值计算出exam2的值，因为在决策边界上的值exam1和exam2存在\(X{{\theta }^{T}}=0\)的关系，因此，我们可以将关系写成\(x_3 = -(\theta_1 + \theta_2 x_2) / \theta_3\)，进而可求解出exam2所对应的值，从而绘制分界线。
此外，我们还需要在这张图中画出原数据的散点图，方便我们对此进行观察。
# 绘制分界线
fig2, ax2 = plt.subplots(figsize=(6, 5))
x = np.linspace(data['exam1'].min(), data['exam1'].max(), 100)
xx, yy = np.meshgrid(x, x)
df = pd.DataFrame()
df['exam1'] = xx.ravel()
df['exam2'] = yy.ravel()
df.insert(0, 'Ones', 1)
df = df.values
# theta1+ theta2 x2 + theta3 x3 = 0 => x3 = -(theta1 + theta2 x2) / theta3
#y_p = -(theta[0][0] + theta[0][1] * x) / theta[0][2]
z = df @ theta.T
print(z.shape)
z = z.reshape(xx.shape)
print(z)
ax2.scatter(positive['exam1'], positive['exam2'], c='r', label='positive')
ax2.scatter(negative['exam1'], negative['exam2'], c='black', label='negative')
# ax2.plot(x, y_p, c='b', label='border')
plt.contour(xx, yy, z, 0)
plt.xlabel("exam1 score")
plt.ylabel("exam2 score")
plt.legend()
从这张图，我们也可以看出总共有100个数据点，只有11个点没有被正确分类，因此，该逻辑分类的训练准确度为89%。
练习2 正则化逻辑回归
在训练的第二部分，我们将要通过加入正则项提升逻辑回归算法。如果你对正则化有点眼生，或者喜欢这一节的方程的背景，请参考在"exercises"文件夹中的"ex2.pdf"。简而言之，正则化是成本函数中的一个术语，它使算法更倾向于“更简单”的模型（在这种情况下，模型将更小的系数）。这个理论助于减少过拟合，提高模型的泛化能力。
在训练的第一部分时，通过观察数据，我们可以看到对此数据的分类模型，分界线大概为直线，而对于看起来非线性的分类模型，我们需要构造多项式特征进行数据训练，此外，还需要尽可能地减少系数，从而避免过拟合，我们可以在损失函数中加入正则化。
接下来，将进行练习2，其题目为：设想你是工厂的生产主管，你有一些芯片在两次测试中的测试结果。对于这两次测试，你想决定是否芯片要被接受或抛弃。为了帮助你做出艰难的决定，你拥有过去芯片的测试数据集，从其中你可以构建一个逻辑回归模型。
数据可视化
和练习1类似的，我们首先也需要进行数据读取之后进行数据可视化
# ex2
# 读取数据
path = "D:\gitcode\Coursera-ML-AndrewNg-Notes-master\Coursera-ML-AndrewNg-Notes-master\code\ex2-logistic regression\ex2data2.txt"
data = pd.read_csv(path, names=["test1", "test2", "accepted"])
print(data.head())
test1 test2 accepted
0 0.051267 0.69956 1
1 -0.092742 0.68494 1
2 -0.213710 0.69225 1
3 -0.375000 0.50219 1
4 -0.513250 0.46564 1
# 可视化数据
positive = data[data['accepted'].isin([1])]
negative = data[data['accepted'].isin([0])]
fig, ax = plt.subplots(figsize=(6, 5))
ax.scatter(positive['test1'], positive['test2'], c='b', label='admitted')
ax.scatter(negative['test1'], negative['test2'], c='r', label='not admitted')
plt.xlabel("test1")
plt.ylabel("test2")
plt.legend()
plt.show()
可以看到这个数据看起来可比前一次的复杂得多。特别地，你会注意到其中没有线性决策界限，来良好的分开两类数据。一个方法是用像逻辑回归这样的线性技术来构造从原始特征的多项式中得到的特征。让我们通过创建一组多项式特征入手吧。
多项式特征转化
在这里我们需要解释一下是如何进行构造的，首先，需要了解是如何通过将特征进行转化成多项式特征的，如果，你要转化的维度是5，那么他构造多项式次数为0-5的式子，例如，如果原本特征为2，分别为\(x_1\)，\(x_2\)，那么，
对于多项式次数为0只能构造\(1\)这个式子，
对于多项式次数为1可以构造\(x_1、x_2\)两个式子，
对于多项式次数为2可以构造\(x_1^2、x_1x_2、x_2^2\)这三个式子，同理类推。
在这里，我们可以自己写特征多项式转化的式子
# 特征映射多项式转化
def feature_mapping(x, y, power):
'''
:param x: 特征1 numpy数组
:param y: 特征2 numpy数组
:param power: 次数
:return: dataframe多项式特征
'''
df = pd.DataFrame()
for p in range(0, power + 1):
for i in range(0, p + 1):
df["f{}{}".format(i, p - i)] = np.power(x, i) * np.power(y, p - i)
return df
初始化数据
在这里，我们将数据进行多项式转化，并进行结果的输出，我们在这里先进行一个变量初始化过程，将数都转化成numpy数组，同时创建theta向量，在前面体到了tnc所要求x0为向量，因此theta需要创建为向量形式。经过多项式转化的式子不会存在全是1的列，因此，我们不在进行插入全为1的列。
# 初始化变量
X = feature_mapping(data["test1"], data["test2"], 6)
print(X.head())
X = X.values
y = data["accepted"].values.reshape(-1, 1)
theta = np.zeros(X.shape[1])
f00 f01 f10 ... f42 f51 f60
0 1.0 0.69956 0.051267 ... 0.000003 2.477505e-07 1.815630e-08
1 1.0 0.68494 -0.092742 ... 0.000035 -4.699318e-06 6.362953e-07
2 1.0 0.69225 -0.213710 ... 0.001000 -3.085938e-04 9.526844e-05
3 1.0 0.50219 -0.375000 ... 0.004987 -3.724126e-03 2.780914e-03
4 1.0 0.46564 -0.513250 ... 0.015046 -1.658422e-02 1.827990e-02
[5 rows x 28 columns]
可以看到转化后具有28个特征。
当然，我们也可以利用sklearn中封装好的函数PolynomialFeatures进行多项式特征转化，同样的我们设置需要转化成最高次数为6，但是对于这个函数返回的为numpy数组和我们所定义的函数的返回值不一样。
from sklearn.preprocessing import PolynomialFeatures
#
poly = PolynomialFeatures(degree=6)
X = poly.fit_transform(data.drop(columns="accepted"))
正则化代价函数
因为进行了多项式特征转化，可能会出现为了防止出现过拟合现象，我们需要使用正则化，因此代价函数，需要进行修改成正则化代价函数。
regularized cost（正则化代价函数）为以下：
def costReg(theta, X, y, lamda):
'''
:param theta: n维向量
:param X: numpy数组 m行n列
:param y: numpy数组 m行1列
:param lamda: 数值 代表正则化学习率
:return: 数值
'''
theta = np.reshape(theta, (1, theta.shape[0]), order='C')
theta1 = theta[:, 1:]
m = X.shape[0]
first = (-y) * np.log(sigmoid(X @ theta.T))
second = (1 - y) * np.log(1 - sigmoid(X @ theta.T))
third = lamda / 2 * (theta1 @ theta1.T)
return (np.sum(first - second) + third[0][0]) / m
在这里，需要请注意等式中的"reg" 项。还注意到另外的一个“学习率”参数。这是一种超参数，用来控制正则化项。我们测试以下，正则化代价函数是否可以使用，其得到最初损失值约为0.69。
print(costReg(theta, X, y, 1))
0.6931471805599454
正则化梯度
同样的，如果我们要使用梯度下降法令这个代价函数最小化，因为我们未对\({{\theta }_{0}}\) 进行正则化，所以梯度下降算法将分两种情形：
\(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]+\frac{\lambda}{2 m} \sum_{j=1}^{n} \theta_{j}^{2}\)
对上面的算法中 j=1,2,...,n 时的更新式子进行调整可得：
\({{\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)}}\)
在这里，我们同样只进行计算梯度即可，在后续使用tnc进行最优化求解
# 计算梯度
def gradientReg(theta, X, y, lamda):
'''
:param theta: n维向量
:param X: numpy数组
:param y: numpy数组
:param lamda:正则化学习率
:return: n维梯度向量
'''
theta = np.reshape(theta, (1, theta.shape[0]), order='C')
theta1 = theta
theta1[0][0] = 0
m = X.shape[0]
return (((sigmoid(X @ theta.T) - y).T @ X + (theta1 * lamda)) / m).flatten()
同样的，我们测试一下最开始的梯度维多少。
print(gradientReg(theta, X, y, 1))
[8.47457627e-03 7.77711864e-05 1.87880932e-02 3.76648474e-02
1.15013308e-02 5.03446395e-02 2.34764889e-02 8.19244468e-03
7.32393391e-03 1.83559872e-02 3.93028171e-02 3.09593720e-03
1.28600503e-02 2.23923907e-03 3.93486234e-02 3.10079849e-02
4.47629067e-03 5.83822078e-03 3.38643902e-03 4.32983232e-03
1.99707467e-02 3.87936363e-02 1.37646175e-03 7.26504316e-03
4.08503006e-04 6.31570797e-03 1.09740238e-03 3.10312442e-02]
SciPy's truncated newton（TNC）实现寻找最优参数
接下来和练习1类似的使用tnc寻找最优参数
# result = opt.fmin_tnc(func=computeCost, x0=theta, fprime=gradient, args=(X, y))
res = opt.fmin_tnc(func=costReg, x0=theta, fprime=gradientReg, args=(X, y, 1))
theta = res[0]
print(res)
(array([ 0. , 0.39178191, -0.00602677, -0.14868226, -0.16345795,
-0.45600957, -0.03803827, -0.07950307, -0.05165685, -0.08286211,
-0.27166958, -0.04280454, -0.12467639, -0.02030656, -0.39959039,
-0.18084981, -0.05335968, -0.05234922, -0.02632121, -0.04113253,
-0.14715525, -0.29410667, -0.01909381, -0.07169553, -0.00138271,
-0.0674712 , -0.00485562, -0.31628257]), 96, 1)
计算精确度
接下来，我们使用练习1所写的预测函数进行精度计算
prediction = predict(theta, X)
num = 0
for i in range(0, len(prediction)):
if (prediction[i] == y[i][0]):
num = num + 1
acc = num / len(X)
print(acc)
0.6610169491525424
可以看到，精确度约为66%
绘制决策边界
在这里，我们一样进行绘制决策边界，其也是寻找等高线为0的线。
fig2, ax2 = plt.subplots(figsize=(6, 5))
x = np.linspace(-1, 1.5, 500)
xx, yy = np.meshgrid(x, x)
z = feature_mapping(xx.ravel(), yy.ravel(), 6).values
z = z @ theta
z = z.reshape(xx.shape)
ax2.scatter(positive['test1'], positive['test2'], c='b', label='admitted')
ax2.scatter(negative['test1'], negative['test2'], c='r', label='not admitted')
plt.legend()
plt.contour(xx, yy, z, 0)
plt.ylim(-.8, 1.2)
plt.show()
可以看到此效果确实没有特别出色
sklearn方法实现正则化逻辑回归
我们也可以使用sklearn提供的方法进行练习2逻辑回归的实现，我们也需要提前利用PolynomialFeatures进行多特征转化，之后利用LogisticRegression提供的l2正则化方法，以及设置正则化学习率为1.0进行模型的训练，最后绘制边界，在这里需要注意的是因为训练完的模型model.coef_为参数矩阵，而model.intercept_为偏置系数，我们需要寻找的等高线应该为-model.intercept_那条。
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import PolynomialFeatures
poly = PolynomialFeatures(degree=6)
X = poly.fit_transform(data.drop(columns="accepted"))
# # X = data.drop(columns=['accepted'])
y = data['accepted']
model = LogisticRegression(penalty='l2', C=1.0)
model.fit(X, y.ravel())
score = model.score(X, y)
theta = model.coef_
print(score)
# 绘制决策边界
x = np.linspace(-1, 1.5, 500)
xx, yy = np.meshgrid(x, x)
df = pd.DataFrame()
df["test1"] = xx.ravel()
df["test2"] = yy.ravel()
z = feature_mapping(xx.ravel(), yy.ravel(), 6).values
fig2, ax3 = plt.subplots(figsize=(6, 5))
z = z @ theta.T
z = z.reshape(xx.shape)
ax3.scatter(positive['test1'], positive['test2'], c='b', label='admitted')
ax3.scatter(negative['test1'], negative['test2'], c='r', label='not admitted')
plt.legend()
plt.contour(xx, yy, z, -model.intercept_)
plt.show()
0.8305084745762712
我们得到模型成绩为0.83是较为好的，通过绘制分界线，也可以通过sklearn训练出来的确实比我们自己训练出来的更优异。
浙公网安备 33010602011771号