补充正则化、绘图以及特征映射内容。
# 序
上则贴子举的例子不含高次项特征，因此绘出的决策边界仅为一条直线，故未讲述如何绘制数据与决策边界。本贴将补充上述未涵盖的内容，并实现 logistic regression 的正则化。
# import modules
import numpy
import pandas
from matplotlib import pyplot# map feature
想要更好地拟合数据，其中一种方法是从每个数据点中创造更多特征。在 map_feature() 函数中将把特征映射到六次方多项式。 stack() 、 concatenate() 与 h/vstack() 的示例将在 Numpy学习记录 中补充。
def map_feature(X_map, degree=6):
"""
:param X_map: (ndarray Shape(m, n)) dataset, m samples by n features
:param degree: the degree of the polynomial you would like to map to
:return: a mapped dataset with more features of higher degree
"""
out = []
for i in range(1, degree+1):
for j in range(i+1):
out.append((X_map[:, 0]**(i-j) * (X_map[:, 1]**j)))
return numpy.stack(out, axis=1)# sigmoid function
后续计算损失、梯度时会经常使用该激活函数，因此不妨构造该函数以便后续调用。激活函数的公式如下：
代码实现如下所示。需要注意的是，参数 z 需要传入 wx+b 。
def sigmoid(z):
return 1 / (1 + numpy.exp(-z))# compute regularized cost
在原始损失函数末尾加上正则化项便得到正则化 logistic regression 的损失计算公式：
代码实现如下：
def compute_cost_reg(X, y, w_in, b_in, lambda_=0.):
"""
compute the cost over all examples
:param X: (ndarray Shape(m, n)) dataset, m samples by n features
:param y: (array_like Shape(m,)) target values for all samples
:param w_in: (array_like Shape(n,)) values of parameters of the model
:param b_in: scalar Values of bias parameter of the model
:param lambda_: the parameter of regularization
:return: the regularized loss
"""
m, n = X.shape
cost_without_reg = 0
for i in range(m):
cost_without_reg += (y[i] * numpy.log(sigmoid(numpy.dot(w_in, X[i]) + b_in)) + (1 - y[i]) * numpy.log(1 - sigmoid(numpy.dot(w_in, X[i]) + b_in)))
cost_without_reg = (-1 / m) * cost_without_reg
reg_cost = 0
for j in range(n):
reg_cost += w_in[j]**2
reg_cost = (lambda_ / (2 * m)) * reg_cost
return cost_without_reg + reg_cost# compute regularized gradient
在梯度下降过程中经历数次迭代，其中每次迭代需重新计算梯度。不妨构造一枚计算梯度的函数以便调用并简化梯度下降的代码。
其中，计算正则化梯度的数学公式如下：
代码实现如下：
def compute_gradient(X, y, w_in, b_in, lambda_=0):
"""
compute gradient for each iteration
:param X: (ndarray Shape(m, n)) dataset, m samples by n features
:param y: (array_like Shape(m,)) target values for all samples
:param w_in: (array_like Shape(n,)) values of parameters of the model
:param b_in: scalar Values of bias parameter of the model
:param lambda_: the parameter of regularization
:return: the gradient dj_dw and dj_db
"""
m, n = X.shape
dj_dw = numpy.zeros(n)
dj_db = 0
# 计算w的梯度
for j in range(n):
for i in range(m):
dj_dw[j] += (sigmoid(numpy.dot(w_in, X[i]) + b_in) - y[i]) * X[i][j]
dj_dw /= m
# 为每一个特征的梯度末尾添加正则化项的偏导数
for j in range(n):
dj_dw[j] += (lambda_ / m) * w_in[j]
# 计算截距项的梯度
for i in range(m):
dj_db += (sigmoid(numpy.dot(w_in, X[i]) + b_in) - y[i])
dj_db /= m
return dj_dw, dj_db# gradient descent
传入映射后的数据集、模型初始参数、学习率、迭代次数及正则化系数并调用正则化损失计算函数、激活函数与正则化梯度计算函数进行梯度下降，最后返回学习得到的模型参数及最后的损失。其中，梯度下降同步更新数学公式如下，较非正则化回归多了正则化项：
代码实现如下：
def gradient_descent(X, y, w_in, b_in, alpha, iters, lambda_=0):
"""
run gradient descent
:param X: (ndarray Shape(m, n)) dataset, m samples by n features
:param y: (array_like Shape(m,)) target values for all samples
:param w_in: (array_like Shape(n,)) values of parameters of the model
:param b_in: scalar value of bias parameter of the model
:param alpha: learning rate α
:param iters: iterations for gradient descent
:param lambda_: the parameter of regularization
:return:
"""
m_in, n_in = X.shape
for inum in range(iters):
dj_dw, dj_db = compute_regularized_gradient(X, y, w_in, b_in, lambda_)
for j in range(n_in):
w_in[j] = w_in[j] - alpha * dj_dw[j]
b_in = b_in - alpha * dj_db
loss_in = compute_regularized_cost(X, y, w_in, b_in, lambda_)
return w_in, b_in, loss_in# predict
正则化仅作用于学习模型参数，将映射后的数据集与学习得到的参数传入激活函数中，将激活函数输出结果按如下规则进行分类，得到预测的目标变量。
实现代码如下：
def predict(X_pred, w_pred, b_pred):
"""
make predictions with learned w and b
:param X_pred: data set with m samples and n features
:param w_pred: values of parameters of the model
:param b_pred: scalar value of bias parameter of the model
:return:
"""
predictions = sigmoid(numpy.dot(X_pred, w_pred) + b_pred)
p = [1 if item >= 0.5 else 0 for item in predictions]
return numpy.array(p)# plot data
将正例、负例数据展示到坐标轴上，并以不同符号颜色进行区分。在此运用Boolean list筛选数组内的元素，需要注意的是被筛选数组的维度。数组的筛选方法具体见贴子 Numpy学习记录 。
def plot_data(X_train, y, ax, positive_label='y=1', negative_label='y=0'):
# 正例列表，当y=1时为True，否则为False。
# 需要注意的是shape为(m,1)，用它过滤shape不为(m,1)或(m,)的ndarray时会匹配不上。
# 而这里X_train为(100,2)，因此需要将Boolean list通过.reshape(-1)转为(m,)
positive = y == 1
# 负例列表，当y=0时为True，否则为False。需要注意的是shape为(m,1)
negative = y == 0
# 通过上述Boolean list绘制正例。用pyplot.scatter()也行
ax.plot(X_train[positive.reshape(-1), 0], X_train[positive.reshape(-1), 1], 'k+', label=positive_label)
# 通过上述Boolean list绘制负例
ax.plot(X_train[negative.reshape(-1), 0], X_train[negative.reshape(-1), 1], 'yo', label=negative_label)# plot decision boundary
在此使用 pyplot.contour() 绘制决策边界，函数 meshgrid() 常与等高线函数一同出现，后者旨在生成网格点，为属性值提供坐标。
meshgrid() 的作用如下图所示：
接下来简单介绍 pyplot.contour() ，前两个参数 X=, Y= 决定位置，既可以传X=XX, Y=YY也可以传X=x, Y=y。第三个参数 Z= 决定位置上的值，需要传维度为 (len(y), len(x)) 的属性值数组。 levels= 传入一个包含想要显示为等高线的值的列表，比如这里决策边界的值为0.5，那么可以传入[0.5]。
def plot_decision_boundary(w_db, b_db, ax):
# 画等高线
x1_ax_contour = numpy.linspace(-1, 1.5, 100)
x2_ax_contour = numpy.linspace(-1, 1.5, 100)
values_contour = numpy.zeros((len(x1_ax_contour), len(x2_ax_contour)))
for i in range(len(x1_ax_contour)):
for j in range(len(x2_ax_contour)):
values_contour[i][j] = sigmoid(numpy.dot(map_feature(numpy.array([x1_ax_contour[i], x2_ax_contour[j]]).reshape((1, 2))), w_db) + b_db)
ax.contour([x1_ax_contour, x2_ax_contour], values_contour.T, levels=[0.5], colors='red')最后讲一下values_contour为什么要转置。由于等高线函数传入决定位置的参数为两个一维列表，第一个参数决定X轴上的位置，第二个参数决定Y轴上的位置。若第一个参数的维度为(n,)，第二个参数的维度为(m,)，那么等高线函数认为的网格点如下：
网格点维度为(m, n)，而通过两个循环计算得到的属性值数组的维度为(n, m)。因此需要将后者转置，才能使每一个网格点对应一个正确的属性值。此外， numpy.meshgrid() 函数返回的每一个矩阵仅对应一个坐标轴上的位置。当把返回的所有矩阵组合在一起，每一个表示一个坐标轴上的位置时，才能形成该多维坐标系上确切的点坐标。
# main
从txt文件中导入数据并对数据集特征映射。初始化所有模型参数后，将它们与学习率、迭代次数、正则化系数一同传入梯度下降函数。值得一提的是，正则化回归的梯度下降函数与非正则化的别无二致，区别仅存在于正则化损失计算与正则化梯度计算。
# 导入数据到一个DataFrame
data = pandas.read_csv('data/data_regular.txt', header=None)
# 获取列数（包括输出变量）
colNum = data.shape[1]
# 定义数据集（不含输出变量）
X_train = data.iloc[:, :colNum - 1]
X_train = numpy.array(X_train.values)
# 特征映射
X_mapped = map_feature(X_train)
# 定义输出空间
y_train = data.iloc[:, colNum - 1: colNum]
y_train = numpy.array(y_train.values)
# 获取样本量和特征数
m, n = X_mapped.shape
# 初始化w, b
numpy.random.seed(1)
w_init = numpy.random.rand(X_mapped.shape[1]) - 0.5
b_init = 1.
# 正则化梯度下降求解模型参数
w, b, loss = gradient_descent(X_mapped, y_train, w_init, b_init, alpha=0.01, iters=10000, lambda_=0.01)
print(f'w = {w}\nb = {b}\nloss = {loss}')
# 根据学习得到的参数进行目标变量的预测
pred = predict(X_mapped, w, b)
# 将预测的目标变量与真实目标变量相比，得到准确率
accuracy = numpy.mean(pred == y_train.reshape(-1)) * 100
print(f'Train Accuracy: {accuracy}')
# 初始化figure和axis
fig, ax = pyplot.subplots()
# 展示数据
plot_data(X_train=X_train, y=y_train, ax=ax)
# 绘决策边界
plot_decision_boundary(w, b, ax)
pyplot.show()# 补充
数组的堆并和连接将在 Numpy学习记录 一贴中展示。
完整代码