线性回归代码实现

线性回归是比较常用的模型。本文会简单介绍线性回归的原理，以及如何用代码实现线性回归模型。

什么是线性回归

简单举一个例子。假设我现在需要贷款，想要知道银行会贷多少钱给我，我从一些渠道知道银行贷款额度跟个人年龄和工资相关。例如我获得了一些数据：

工资	年龄	额度
4000	25	20000
5000	28	35000
7500	33	50000

我从这些数据中得到三者之间的关系，建立了一个函数模型，例如： $y = \theta_1x_1 + \theta_2x_2 + \varepsilon$ ，然后我用这个函数模型，输入我自己的工资和年龄，从而预测我可以从银行贷多少钱。这就是线性回归分析，因为包含两个自变量（工资和年龄），所以也叫多元线性回归。其中， $\theta_1$ 和 $\theta_2$ 两个参数表示工资和年龄对贷款额度的影响程度， $\varepsilon$ 是误差项，服从均值为0的正态分布。

一些数学式子

假设有个线性回归问题有m个变量，则可以建立函数 $h_\theta(x) = \theta_0 + \theta_1x_1 + \theta_2x_2 + \cdots + \theta_mx_m$ ，其中 $\theta_0$ 为偏置项。式子可以整合为 $h_\theta(x) = \sum_{i = 0}^{n}\theta_ix_i = \theta^Tx$ 。

因为真实值和预测值之间肯定会存在差异的，所以对于每个样本来说， $y^{(i)} = \theta^Tx^{(i)} + \varepsilon^{(i)}$ 。

误差是服从均值为0的正态分布。所以 $p(\varepsilon^{(i)}) = \frac{1}{\sqrt{2\pi}\sigma}exp(-\frac{(\varepsilon^{(i)})^2}{2\sigma^2})$ 。

最小二乘法： $J(\theta) = \frac{1}{2}\sum_{i = 1}^{m}(y^{(i)} - \theta^Tx^{(i)})^2$

批量梯度下降： $\theta_j = \theta_j - \alpha\frac{1}{10}\sum_{k = i}^{i + 9}(h_\theta(x^{(k)}) - y^{(k)})x_j^{(k)}$ （ $\alpha$ 为学习率，一般较小；每次更新选择10个数据来计算；具体数量根据实际情况调整）

代码实现

数据预处理

写一个prepare_for_training函数，对数据进行函数变换、标准化等操作。最后返回处理过的数据，以及均值和标准差。

def prepare_for_training(data, polynomial_degree=0, sinusoid_degree=0, normalize_data=True):

    # 计算样本总数
    num_examples = data.shape[0]

    data_processed = np.copy(data)

    # 预处理
    features_mean = 0
    features_deviation = 0
    data_normalized = data_processed
    if normalize_data:
        (
            data_normalized,
            features_mean,
            features_deviation
        ) = normalize(data_processed)

        data_processed = data_normalized

    # 特征变换sinusoidal
    if sinusoid_degree > 0:
        sinusoids = generate_sinusoids(data_normalized, sinusoid_degree)
        data_processed = np.concatenate((data_processed, sinusoids), axis=1)

    # 特征变换polynomial
    if polynomial_degree > 0:
        polynomials = generate_polynomials(data_normalized, polynomial_degree)
        data_processed = np.concatenate((data_processed, polynomials), axis=1)

    # 加一列1
    data_processed = np.hstack((np.ones((num_examples, 1)), data_processed))

    return data_processed, features_mean, features_deviation

线性回归模块

写一个LinearRegression类，包含线性回归相关的方法。

    def train(self, alpha, num_iterations = 500):
        """
                    训练模块，执行梯度下降
        """
        cost_history = self.gradient_descent(alpha, num_iterations)
        return self.theta, cost_history

    def gradient_descent(self, alpha, num_iterations):
        """
                    实际迭代模块，会迭代num_iterations次
        """
        cost_history = []
        for _ in range(num_iterations):
            self.gradient_step(alpha)
            cost_history.append(self.cost_function(self.data, self.labels))
        return cost_history

    def gradient_step(self, alpha):
        """
                    梯度下降参数更新计算方法，注意是矩阵运算
        """
        num_examples = self.data.shape[0]
        prediction = LinearRegression.hypothesis(self.data, self.theta)
        delta = prediction - self.labels
        theta = self.theta
        theta = theta - alpha * (1/num_examples)*(np.dot(delta.T, self.data)).T
        self.theta = theta

    def cost_function(self, data, labels):
        """
                    损失计算方法
        """
        num_examples = data.shape[0]
        delta = LinearRegression.hypothesis(data, self.theta) - labels
        cost = np.dot(delta, delta.T)
        return cost[0][0]

    @staticmethod
    def hypothesis(data, theta):
        predictions = np.dot(data, theta)
        return predictions

    def get_cost(self, data, labels):
        data_processed = prepare_for_training.prepare_for_training(data,
                                    self.polynomial_degree,
                                    self.sinusoid_degree,
                                    self.normalize_data)[0]
        return self.cost_function(data_processed, labels)

    def predict(self, data):
        """
                    用训练的参数模型，预测得到回归值结果
        """
        data_processed = prepare_for_training.prepare_for_training(data,
                                    self.polynomial_degree,
                                    self.sinusoid_degree,
                                    self.normalize_data)[0]
        predictions = LinearRegression.hypothesis(data_processed, self.theta)

        return predictions

实际应用

用LinearRegression类进行建模、预测、计算损失等。

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

from linear_regression import LinearRegression

data = pd.read_csv('../data/world-happiness-report-2017.csv')

# 得到训练和测试数据
train_data = data.sample(frac= 0.8)
test_data = data.drop(train_data.index)

input_param_name = 'Economy..GDP.per.Capita.'
output_param_name = 'Happiness.Score'

x_train = train_data[[input_param_name]].values
y_train = train_data[[output_param_name]].values

x_test = test_data[[input_param_name]].values
y_test = test_data[[output_param_name]].values

plt.scatter(x_train, y_train, label='Train data')
plt.scatter(x_test, y_test, label='Test data')
plt.xlabel(input_param_name)
plt.ylabel(output_param_name)
plt.title('Happy')
plt.legend()
plt.show()

num_iterations = 500
learning_rate = 0.01

linear_regression = LinearRegression(x_train, y_train)
(theta, cost_history) = linear_regression.train(learning_rate, num_iterations)

print('开始时的损失：', cost_history[0])
print('训练后的损失：', cost_history[-1])

plt.plot(range(num_iterations), cost_history)
plt.xlabel('Iteration')
plt.ylabel('Cost')
plt.title('GD')
plt.show()

predictions_num = 100
x_predictions = np.linspace(x_train.min(), x_train.max(), predictions_num).reshape(predictions_num,1)
y_predictions = linear_regression.predict(x_predictions)

plt.scatter(x_train, y_train, label='Train data')
plt.scatter(x_test, y_test, label='Test data')
plt.plot(x_predictions, y_predictions, 'r', label='Prediction')
plt.xlabel(input_param_name)
plt.ylabel(output_param_name)
plt.title('Happy')
plt.legend()
plt.show()

运行结果如下。

开始时的损失： 30.901003728555025
训练后的损失： 0.3040912375218431