2021/01/23 19:05

Kumud Lakara作者

从零开始，用Python徒手写线性回归

先放下 Scikit-learn，我们来看一看真正的技术。

对于大多数数据科学家而言，线性回归方法是他们进行统计学建模和预测分析任务的起点。这种方法已经存在了 200 多年，并得到了广泛研究，但仍然是一个积极的研究领域。由于良好的可解释性，线性回归在商业数据上的用途十分广泛。当然，在生物数据、工业数据等领域也不乏关于回归分析的应用。

另一方面，Python 已成为数据科学家首选的编程语言，能够应用多种方法利用线性模型拟合大型数据集显得尤为重要。

如果你刚刚迈入机器学习的大门，那么使用 Python 从零开始对整个线性回归算法进行编码是一次很有意义的尝试，让我们来看看怎么做吧。

数据

机器学习问题的第一步是获取数据，没有可以学习的数据就没有机器学习。本文将使用非常常规的线性回归数据集——房价预测数据集。

这是一个包含俄勒冈州波特兰市房价的简单数据集。该数据集中第一列是房屋面积（以平方英尺为单位），第二列是卧室的数量，第三列是房屋价格。该数据集中有多个特征（例如，house_size 和房间数），因此我们将研究多元线性回归，标签 (y) 是我们将要预测的房价。

首先定义用于加载数据集的函数：

def load_data(filename):
    df = pd.read_csv(filename, sep=",", index_col=False)
    df.columns = ["housesize", "rooms", "price"]
    data = np.array(df, dtype=float)
    plot_data(data[:,:2], data[:, -1])
    normalize(data)
    return data[:,:2], data[:, -1]

我们稍后将调用上述函数来加载数据集。此函数返回 x 和 y。

归一化数据

上述代码不仅加载数据，还对数据执行归一化处理并绘制数据点。在查看数据图之前，我们首先了解上述代码中的 normalize(data)。

查看原始数据集后，你会发现第二列数据的值（房间数量）比第一列（即房屋面积）小得多。该模型不会将此数据评估为房间数量或房屋面积，对于模型来说，它们只是一些数字。机器学习模型中某些列（或特征）的数值比其他列高可能会造成不想要的偏差，还可能导致方差和数学均值的不平衡。出于这些原因，也为了简化工作，我们建议先对特征进行缩放或归一化，使其位于同一范围内（例如 [-1,1] 或 [0,1]），这会让训练容易许多。因此我们将使用特征归一化，其数学表达如下：

Z = (x — μ) / σ
μ : mean
σ : standard deviation

其中 z 是归一化特征，x 是非归一化特征。有了归一化公式，我们就可以为归一化创建一个函数：

def normalize(data):
    for i in range(0,data.shape[1]-1):
        data[:,i] = ((data[:,i] - np.mean(data[:,i]))/np.std(data[:, i]))

上述代码遍历每一列，并使用每一列中所有数据元素的均值和标准差对其执行归一化。

绘制数据

在对线性回归模型进行编码之前，我们需要先问「为什么」。

为什么要使用线性回归解决这个问题？这是一个非常有用的问题，在写任何具体代码之前，你都应该非常清楚要使用哪种算法，以及在给定数据集和待解决问题的情况下，这是否真的是最佳选择。

我们可以通过绘制图像来证明对当前数据集使用线性回归有效的原因。为此，我们在上面的 load_data 中调用了 plot_data 函数，现在我们来定义一下 plot_data 函数：

def plot_data(x, y):
    plt.xlabel('house size')
    plt.ylabel('price')
    plt.plot(x[:,0], y, 'bo')
    plt.show()

调用该函数，将生成下图：

房屋面积与房屋价格关系图。

如上图所示，我们可以粗略地拟合一条线。这意味着使用线性近似能够做出较为准确的预测，因此可以采用线性回归。

准备好数据之后就要进行下一步，给算法编写代码。

假设

首先我们需要定义假设函数，稍后我们将使用它来计算代价。对于线性回归，假设是：

但数据集中只有 2 个特征，因此对于当前问题，假设是：

其中 x1 和 x2 是两个特征（即房屋面积和房间数量）。然后编写一个返回该假设的简单 Python 函数：

def h(x,theta):
    return np.matmul(x, theta)

接下来我们来看代价函数。

代价函数

使用代价函数的目的是评估模型质量。

代价函数的等式为：

代价函数的代码如下：

def cost_function(x, y, theta):
    return ((h(x, theta)-y).T@(h(x, theta)-y))/(2*y.shape[0])

到目前为止，我们定义的所有 Python 函数都与上述线性回归的数学意义完全相同。接下来我们需要将代价最小化，这就要用到梯度下降。

梯度下降

梯度下降是一种优化算法，旨在调整参数以最小化代价函数。

梯度下降的主要更新步是：

因此，我们将代价函数的导数乘以学习率（α），然后用参数（θ）的当前值减去它，获得新的更新参数（θ）。

def gradient_descent(x, y, theta, learning_rate=0.1, num_epochs=10):
    m = x.shape[0]
    J_all = []
    
    for _ in range(num_epochs):
        h_x = h(x, theta)
        cost_ = (1/m)*(x.T@(h_x - y))
        theta = theta - (learning_rate)*cost_
        J_all.append(cost_function(x, y, theta))

    return theta, J_all

gradient_descent 函数返回 theta 和 J_all。theta 显然是参数向量，其中包含假设的θs 值，J_all 是一个列表，包含每个 epoch 后的代价函数。J_all 变量并非必不可少，但它有助于更好地分析模型。

整合到一起

接下来要做的就是以正确的顺序调用函数：

x,y = load_data("house_price_data.txt")
y = np.reshape(y, (46,1))
x = np.hstack((np.ones((x.shape[0],1)), x))
theta = np.zeros((x.shape[1], 1))
learning_rate = 0.1
num_epochs = 50
theta, J_all = gradient_descent(x, y, theta, learning_rate, num_epochs)
J = cost_function(x, y, theta)
print("Cost: ", J)
print("Parameters: ", theta)

#for testing and plotting cost 
n_epochs = []
jplot = []
count = 0
for i in J_all:
    jplot.append(i[0][0])
    n_epochs.append(count)
    count += 1
jplot = np.array(jplot)
n_epochs = np.array(n_epochs)
plot_cost(jplot, n_epochs)

test(theta, [1600, 2])

首先调用 load_data 函数载入 x 和 y 值。x 值包含训练样本，y 值包含标签（在这里就是房屋的价格）。

你肯定注意到了，在整个代码中，我们一直使用矩阵乘法的方式来表达所需。例如为了得到假设，我们必须将每个参数（θ）与每个特征向量（x）相乘。我们可以使用 for 循环，遍历每个样本，每次都执行一次乘法，但如果训练的样本过多，这可能不是最高效的方法。

在这里更有效的方式是使用矩阵乘法。本文所用的数据集具备两个特征：房屋面积和房间数，即我们有（2+1）三个参数。将假设看作图形意义上的一条线，用这种方式来思考额外参数θ0，最终额外的θ0 也要使这条线符合要求。

有利的假设函数图示。

现在我们有了三个参数和两个特征。这意味着θ或参数向量（1 维矩阵）的维数是 (3,1)，但特征向量的维度是 (46,2)。你肯定会注意到将这样两个矩阵相乘在数学上是不可能的。再看一遍我们的假设：

如果你仔细观察的话，实际上这很直观：如果在特征向量 (x) {维度为 (46, 3)} 的开头添加额外的一列，并且对 x 和 theta 执行矩阵乘法，将得出 hθ(x) 的方程。

记住，在实际运行代码来实现此功能时，不会像 hθ(x) 那样返回表达式，而是返回该表达式求得的数学值。在上面的代码中，x = np.hstack((np.ones((x.shape[0],1)), x)) 这一行在 x 开头加入了额外一列，以备矩阵乘法需要。

在这之后，我们用零初始化 theta 向量，当然你也可以用一些小随机值来进行初始化。我们还指定了训练学习率和 epoch 数。

定义完所有超参数之后，我们就可以调用梯度下降函数，以返回所有代价函数的历史记录以及参数 theta 的最终向量。在这里 theta 向量定义了最终的假设。你可能注意到，由梯度下降函数返回的 theta 向量的维度为 (3,1)。

还记得函数的假设吗？

所以我们需要三个θ，theta 向量的维度为 (3,1)，因此 theta [0]、theta [1] 和 theta [2] 实际上分别为θ0、θ1 和 θ2。J_all 变量是所有代价函数的历史记录。你可以打印出 J_all 数组，来查看代价函数在梯度下降的每个 epoch 中逐渐减小的过程。

代价和 epoch 数量的关系图。

我们可以通过定义和调用 plot_cost 函数来绘制此图，如下所示：

def plot_cost(J_all, num_epochs):
    plt.xlabel('Epochs')
    plt.ylabel('Cost')
    plt.plot(num_epochs, J_all, 'm', linewidth = "5")
    plt.show()

现在我们可以使用这些参数来找到标签，例如给定房屋面积和房间数量时的房屋价格。

测试

现在你可以测试调用测试函数的代码，该函数会将房屋面积、房间数量和 logistic 回归模型返回的最终 theta 向量作为输入，并输出房屋价格。

def test(theta, x):
    x[0] = (x[0] - mu[0])/std[0]
    x[1] = (x[1] - mu[1])/std[1]

    y = theta[0] + theta[1]*x[0] + theta[2]*x[1]
    print("Price of house: ", y)

完整代码

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

#variables to store mean and standard deviation for each feature
mu = []
std = []

def load_data(filename):
    df = pd.read_csv(filename, sep=",", index_col=False)
    df.columns = ["housesize", "rooms", "price"]
    data = np.array(df, dtype=float)
    plot_data(data[:,:2], data[:, -1])
    normalize(data)
    return data[:,:2], data[:, -1]

def plot_data(x, y):
    plt.xlabel('house size')
    plt.ylabel('price')
    plt.plot(x[:,0], y, 'bo')
    plt.show()

def normalize(data):
    for i in range(0,data.shape[1]-1):
        data[:,i] = ((data[:,i] - np.mean(data[:,i]))/np.std(data[:, i]))
        mu.append(np.mean(data[:,i]))
        std.append(np.std(data[:, i]))


def h(x,theta):
    return np.matmul(x, theta)

def cost_function(x, y, theta):
    return ((h(x, theta)-y).T@(h(x, theta)-y))/(2*y.shape[0])

def gradient_descent(x, y, theta, learning_rate=0.1, num_epochs=10):
    m = x.shape[0]
    J_all = []
    
    for _ in range(num_epochs):
        h_x = h(x, theta)
        cost_ = (1/m)*(x.T@(h_x - y))
        theta = theta - (learning_rate)*cost_
        J_all.append(cost_function(x, y, theta))

    return theta, J_all 

def plot_cost(J_all, num_epochs):
    plt.xlabel('Epochs')
    plt.ylabel('Cost')
    plt.plot(num_epochs, J_all, 'm', linewidth = "5")
    plt.show()

def test(theta, x):
    x[0] = (x[0] - mu[0])/std[0]
    x[1] = (x[1] - mu[1])/std[1]

    y = theta[0] + theta[1]*x[0] + theta[2]*x[1]
    print("Price of house: ", y)

x,y = load_data("house_price_data.txt")
y = np.reshape(y, (46,1))
x = np.hstack((np.ones((x.shape[0],1)), x))
theta = np.zeros((x.shape[1], 1))
learning_rate = 0.1
num_epochs = 50
theta, J_all = gradient_descent(x, y, theta, learning_rate, num_epochs)
J = cost_function(x, y, theta)
print("Cost: ", J)
print("Parameters: ", theta)

#for testing and plotting cost 
n_epochs = []
jplot = []
count = 0
for i in J_all:
    jplot.append(i[0][0])
    n_epochs.append(count)
    count += 1
jplot = np.array(jplot)
n_epochs = np.array(n_epochs)
plot_cost(jplot, n_epochs)

test(theta, [1600, 3])

总结

这就是线性回归的全部代码了。

现在你已经学会了从零开始成功编写线性回归模型。能够理解和编写整个算法并不是一件容易的事，你或许需要时不时地回看才能完全理解。但这些努力是值得的，线性回归通常是人们学习机器学习算法的第一步，在这之后你可以选择另一个适用于线性回归处理的数据集，并尝试刚写好的算法。

原文链接：

https://towardsdatascience.com/coding-linear-regression-from-scratch-c42ec079902

理论线性回归Python

相关数据

机器学习技术

机器学习是人工智能的一个分支，是一门多领域交叉学科，涉及概率论、统计学、逼近论、凸分析、计算复杂性理论等多门学科。机器学习理论主要是设计和分析一些让计算机可以自动“学习”的算法。因为学习算法中涉及了大量的统计学理论，机器学习与推断统计学联系尤为密切，也被称为统计学习理论。算法设计方面，机器学习理论关注可以实现的，行之有效的学习算法。

来源：Mitchell, T. (1997). Machine Learning. McGraw Hill.

参数技术

在数学和统计学裡，参数（英语：parameter）是使用通用变量来建立函数和变量之间关系（当这种关系很难用方程来阐述时）的一个数量。

来源：维基百科

数据科学技术

数据科学，又称资料科学，是一门利用数据学习知识的学科，其目标是通过从数据中提取出有价值的部分来生产数据产品。它结合了诸多领域中的理论和技术，包括应用数学、统计、模式识别、机器学习、数据可视化、数据仓库以及高性能计算。数据科学通过运用各种相关的数据来帮助非专业人士理解问题。

来源：维基百科

学习率技术

在使用不同优化器（例如随机梯度下降，Adam）神经网络相关训练中，学习速率作为一个超参数控制了权重更新的幅度，以及训练的速度和精度。学习速率太大容易导致目标（代价）函数波动较大从而难以找到最优，而弱学习速率设置太小，则会导致收敛过慢耗时太长

来源：Liu, T. Y. (2009). Learning to rank for information retrieval. Foundations and Trends® in Information Retrieval, 3(3), 225-331. Wikipedia

超参数技术

在机器学习中，超参数是在学习过程开始之前设置其值的参数。相反，其他参数的值是通过训练得出的。不同的模型训练算法需要不同的超参数，一些简单的算法（如普通最小二乘回归）不需要。给定这些超参数，训练算法从数据中学习参数。相同种类的机器学习模型可能需要不同的超参数来适应不同的数据模式，并且必须对其进行调整以便模型能够最优地解决机器学习问题。在实际应用中一般需要对超参数进行优化，以找到一个超参数元组（tuple），由这些超参数元组形成一个最优化模型，该模型可以将在给定的独立数据上预定义的损失函数最小化。

来源：Wikipedia

导数技术

导数（Derivative）是微积分中的重要基础概念。当函数y=f(x)的自变量x在一点x_0上产生一个增量Δx时，函数输出值的增量Δy与自变量增量Δx的比值在Δx趋于0时的极限a如果存在，a即为在x0处的导数，记作f'(x_0) 或 df(x_0)/dx。

来源：百度百科

线性回归技术

在现实世界中，存在着大量这样的情况：两个变量例如X和Y有一些依赖关系。由X可以部分地决定Y的值，但这种决定往往不很确切。常常用来说明这种依赖关系的最简单、直观的例子是体重与身高，用Y表示他的体重。众所周知，一般说来，当X大时，Y也倾向于大，但由X不能严格地决定Y。又如，城市生活用电量Y与气温X有很大的关系。在夏天气温很高或冬天气温很低时，由于室内空调、冰箱等家用电器的使用，可能用电就高，相反，在春秋季节气温不高也不低，用电量就可能少。但我们不能由气温X准确地决定用电量Y。类似的例子还很多，变量之间的这种关系称为“相关关系”，回归模型就是研究相关关系的一个有力工具。

来源：王松桂等编线性统计模型线性回归与方差分析高等教育出版社 Wikipedia

梯度下降技术

梯度下降是用于查找函数最小值的一阶迭代优化算法。要使用梯度下降找到函数的局部最小值，可以采用与当前点的函数梯度（或近似梯度）的负值成比例的步骤。如果采取的步骤与梯度的正值成比例，则接近该函数的局部最大值，被称为梯度上升。

来源：Vapnik V. N. (2000). The Nature of Statistical Learning Theory. Information Science and Statistics. Springer-Verlag.Wikipedia

回归分析技术

回归分析是一种用于估计变量之间的关系（当一个自变量变化而其它变量固定时，因变量会如何变化）的统计过程，在预测任务中有广泛的应用。回归分析模型有不同的种类，其中最流行的是线性回归和逻辑回归（Logistic Regression）。另外还有多变量回归、泊松回归、逐步回归、脊回归（Ridge Regression）、套索回归（Lasso Regression）和多项式回归等等。随机梯度下降（SGD）就是一种起源于回归分析的常用方法，可用于控制复杂度。

来源：机器之心