Numpy sigmoid回归

Numpy和Scipy是科学计算的Python库，如今已成为科技工作者的不可或缺的工具。本文将介绍其中的一个主题，即Numpy中的sigmoid回归。

阅读更多：Numpy 教程

Sigmoid回归

回归分析是数据分析领域中的重要方法，它可以用来预测未来的数据，并寻找数据之间的相互关系。一种常见的回归方法是sigmoid回归。Sigmoid函数是一个S形曲线，它在区间(-∞,+∞)之间都有定义。在分类问题中，Sigmoid函数可以将不同的数据集分为两个类别。
Sigmoid是一个反函数，本身可以用于逆转变换，即从一个值域到另一个值域。

Numpy中的Sigmoid函数

Numpy的sigmoid函数可以用如下代码定义，其中x是一个向量，y是其对应的sigmoid值：

import numpy as np
def sigmoid(x):
    return 1 / (1 + np.exp(-x))

y = sigmoid(x)

下面，我们将使用如上的sigmoid方法，来为一些数据进行分类。

假设我们有如下数据，其中前两列是x,y坐标轴的特征，最后一列是目标分类结果：

x	y	target
1.0	2.1	0
2.0	1.4	0
1.1	1.0	1
1.3	1.1	1

我们将会使用numpy对数据进行处理，并通过scipy中的optimize工具包来达到我们的目标——对目标分类结果进行预测。

数据处理

大多数时候，我们的输入数据都需要进行一些简单的处理。在python中，我们可以使用numpy库方便地对输入数据进行矩阵运算。首先，我们需要将数据加载到numpy中，并将目标分类结果进行处理，使其成为1或0的向量：

import numpy as np
import matplotlib.pyplot as plt

data = np.loadtxt('data.dat', delimiter=',')
x = np.array(data[:, 0:2])
y = np.array(data[:, 2])
y[y == 0] = -1  # 针对y==0的情况，将其修改为-1

定义代价函数

代价函数是用于描述模型的预测值与实际值之间误差的函数。在sigmoid回归中，我们采用的是逻辑代价函数，它可以描述模型误差的程度，并帮助我们对模型进行优化。逻辑代价函数的式子如下：
$J(w)=\sum\limits_{i=1}^{m}ln(1+e^{-y^{(i)}w^Tx^{(i)}})$
其中，w是一个向量，x是输入的样本矩阵，y是目标向量。

在Numpy中，我们可以用如下代码定义已知数据集的代价函数：

import numpy as np

def cost(theta, x, y):
    h = sigmoid(np.dot(x, theta))
    return np.sum(np.log(1 + np.exp(- y * h)))

确定优化参数

确定优化参数的过程，就是利用梯度下降来寻找最优解的过程。我们将使用SciPy的优化算法来求解参数。SciPy最常用的优化函数之一便是scipy.optimize.minimize函数。
使用scipy.optimize.minimize函数需要定义一个损失函数，这个损失函数需要有一个输入参数并且根据这个参数来计算损失。在此过程中，我们需要使用到Numpy库的一些线性代数和矩阵运算的函数：

from scipy.optimize import minimize

def fit(x, y, initial_theta):
    res = minimize(fun=cost, x0=initial_theta, args=(x,y), method='TNC', jac=gradient)
    return res.x

# 初始化参数
initial_theta = np.zeros(x.shape[1])
theta = fit(x, y, initial_theta)

预测结果的验证

接下来，我们需要使用得到的参数θ来进行预测结果的验证，并将其绘制出来。具体地，我们将把横坐标轴作为样本的特征点x，纵坐标轴作为样本的特征点y，用红和蓝颜色的点分别表示两个不同分类结果：

# 预测结果的可视化
plt.scatter(x[y == 1, 0], x[y == 1, 1], marker='o', label='class_1')
plt.scatter(x[y == -1, 0], x[y == -1, 1], marker='x', label='class_2')

# 绘制决策边界
plot_x = np.array([np.min(x[:, 0]), np.max(x[:, 0])])
plot_y = (- 1.0 / theta[1]) * (theta[0] * plot_x + theta[2])
plt.plot(plot_x, plot_y)

plt.legend()
plt.show()

执行上述代码，我们可以得到如下的预测结果可视化图：

可以看出，我们的分类模型在两类数据之间建立了合理的边界。