机器学习——线性回归模型及python代码实现

《机器学习：公式推导与代码实践》鲁伟著读书笔记。在机器学习的学习过程中，相信大家首先要学习的就是线性模型。而线性模型中，线性回归（Linear Regression）是一种非常经典方法。现在我从线性回归的数学原理出发，手推数学公式，并结合python代码，对线性回归模型进行系统性的总结。

线性回归的数学原理

线性回归从实质上说为通过训练进而得到一个线性模型来根据输入数据$X$来拟合输出$y$。面对多元化的线性回归问题，设训练模型的输入矩阵为$X\in {\mathbb{R}}^{m\times d}$，其中$m$为训练集输入数据的组数，$d$为变量$X$的特征维数。输出$y\in {\mathbb{R}}^{m\times 1}$为了方便矩阵化的最小二乘法的推导，可以设权重为$\omega \in {\mathbb{R}}^{d\times 1}$与偏置为$b\in {\mathbb{R}}^{m\times 1}$。我采用均方误差来度量预测与标签之间的损失，所以我们的线性回归任务的优化目标便为使线性回归模型的预测值和输出的真实值之间的均方误差最小化。
参数优化的目标函数为：$$({\omega }^{\ast },b^{\ast })=\arg \min (y-(X\omega +b){)}^T(y-(X\omega +b))$$令$L=(y-(X\omega +b){)}^T(y-(X\omega +b))$，并就参数优化的目标函数对$\omega$进行求导，求得权重的梯度。推导过程如下：$$L=y^{T}y-y^T(X\omega +b)-\left({\omega }^T{X}^T+b^T\right)y+{\omega }^T{X}^{T}X\omega +{\omega }^T{X}^{T}b+b^{T}X\omega +b^{T}b$$又因为矩阵微分公式：$\frac{\partial a^{T}x}{\partial x}=\frac{\partial x^{T}a}{\partial x}=a$，$\frac{\partial x^{T}Ax}{\partial x}=(A+A^T)x$。
$$\frac{\partial L}{\partial \omega }=\frac{\partial y^{T}y}{\partial \omega }-\frac{\partial y^T(X\omega +b)}{\partial \omega }-\frac{\partial \left({\omega }^T{X}^T+b^T\right)y}{\partial \omega }+\frac{\partial {\omega }^T{X}^{T}X\omega }{\partial \omega }+\frac{\partial {\omega }^T{X}^{T}b}{\partial \omega }+\frac{\partial b^{T}X\omega }{\partial \omega }+\frac{\partial b^{T}b}{\partial \omega }$$通过化简得到：
$$\frac{\partial L}{\partial \omega }=0-X^{T}y-X^{T}y+(X^{T}X+X^{T}X)\omega +X^{T}b+X^{T}b+0$$$$\frac{\partial L}{\partial \omega }=2X^T(X\omega +b-y)$$即：$$\frac{\partial L}{\partial \omega }=2X^T(\hat{y}-y)$$下一步对偏置进行求导，得到偏置的梯度。$$\frac{\partial L}{\partial b}=\frac{\partial y^{T}y}{\partial b}-\frac{\partial y^T(X\omega +b)}{\partial b}-\frac{\partial \left({\omega }^T{X}^T+b^T\right)y}{\partial b}+\frac{\partial {\omega }^T{X}^{T}X\omega }{\partial b}+\frac{\partial {\omega }^T{X}^{T}b}{\partial b}+\frac{\partial b^{T}X\omega }{\partial b}+\frac{\partial b^{T}b}{\partial b}$$$$\frac{\partial L}{\partial b}=0-y-y+X\omega +X\omega +2b$$$$\frac{\partial L}{\partial b}=2(X\omega +b-y)$$即：$$\frac{\partial L}{\partial b}=2(\hat{y}-y)$$最后便可以运用权重和偏置的梯度值完成在训练过程中参数的更新，例如：${\omega }^{\ast }=\omega -learningrate\ast \frac{\partial L}{\partial \omega }$；$b^{\ast }=b-learningrate\ast \frac{\partial L}{\partial b}$。

线性回归的NumPy手撕代码

按照推理得出的参数更新梯度，线性回归模型主体部分的实现如下：

初始化权重和偏置参数

def init_params(train_dim):
	'''
	输入：
	train_dim：样本特征数
	输出：
	w：初始化后的权重
	b：初始化后的偏置
	'''
	w = np.zeros((train_dim,1))
	b = 0
	return w,b

线性回归算法

def linear_regress(X,y,w,b):
	'''
	输入：
	X：输入数据
	y：输出标签
	w：权重
	b：偏置
	输出：
	y_hat：预测值
	loss：预测值与真实值的均方误差
	dw：权重的一阶梯度
	db：偏置的一阶梯度
	'''
	num_train = X.shape[0]
	num_feature = X.shape[1]
	y_hat = np.dot(X,w) + b
	loss = np.sum((y_hat-y)**2)/num_train
	dw = np.dot(X.T,(y_hat-y))/num_train
	db = np.sum((y_hat-y))/num_train
	return y_hat, loss, dw, db	

定义训练过程

采用梯度下降的方法对参数进行更新。

def train(X, y, learning_rate=0.01, epochs=10000):
	'''
	输入：
	X：输入数据
	y：输出标签
	learning_rate：学习率
	epochs：迭代次数
	输出：
	loss_his：每一代的误差
	params：参数字典
	grads：优化后的梯度
	'''
	loss_his = []
	w, b = init_params(X.shape[1])
	for i in range(epochs):
		y_hat, loss, dw, db = linear_regress(X, y, w, b)
		w += -learning_rate*dw
		b += -learning_rate*db
		loss_his.append(loss)
	params = {'w':w, 'b':b}
	grads = {'dw':dw,'db':db}
	return loss_his, params, grads

定义预测过程

def predict(X, params):
	'''
	输入：
	X：测试数据集
	params：模型训练参数
	输出：
	y_pre：预测值
	'''
	w = params['w']
	b = params['b']
	y_pre = np.dot(X, w) + b
	return y_pre

实战

我们采用sklearn中的diabetes数据集进行测试。其一共有442个样本量，每组样本有10个特征。将数据集按照8：2 的比例进行划分训练集和测试集。具体训练过程如下：
导入数据并划分数据集

from sklearn.datasets import load_diabetes
from sklearn.utils import shuffle # 打乱数据的函数
import numpy as np

diabetes = load_diabetes()
data,target = diabetes.data, diabetes.target
X, y = shuffle(data, target, random_state=13) # 固定random_state后，每次打乱的数据集都是一样的
train_data = int(X.shape[0]*0.8)
X_train, y_train = X[:train_data], y[:train_data]
X_test, y_test = X[train_data:], y[train_data:]
y_train = y_train.reshape((-1,1)) # 将target从(1，m)转变为(m,1)，其实就是转置了一下
y_test = y_test.reshape((-1,1))

# train model
loss_his, params, grads = train(X_train, y_train, 0.01, 10000)

# predict model
y_pre = predict(X_test, params)

# 评价模型
def r2_score(y_test, y_pre):
	y_avg = np.mean(y_test)
	ss_tot = np.sum((y_test-y_avg)**2) # 总离差平方和
	ss_res = np.sum((y_test-y_pre)**2) # 残差平方和
	r2 = 1-(ss_res/ss_tot)
	return r2
pre_r2 = r2_score(y_test, y_pre)
print(pre_r2)

最后计算得到的线性回归模型$R^2$系数为0.5334188778093094。预测精准度不是很高，说明线性回归的预测性能不高。下一个章节，我将总结对数几率回归模型，即logistic regression。