在统计学中，线性回归(Linear Regression)是利用称为线性回归方程的最小平方函数对一个或多个自变量和因变量之间关系进行建模的一种回归分析。其中只有一个自变量的情况称为简单回归,大于一个自变量情况的叫做多元回归。
一、梯度下降算法实现线性回归
首先我们载入奥运马拉松数据。这是从1896年第一届奥运会到伦敦奥运会男子马拉松比赛的奥运马拉松时间数据。
import podsdata = pods.datasets.olympic_marathon_men()x = data['X']y = data['Y']
并打印数据集
importpylabaspltplt.plot(x, y, 'rx')plt.xlabel('year')plt.ylabel('pace in min/km')
现在，我们将采用我们最大似然法，将一条线
与绘制的数据拟合，并尽量减少误差函数：
关于其中变量m和c，我们可以先初始化一个假设值m，之后利用最大似然法迭代更新从而得到c。
m = -0.4c = 80
对于c求偏导，可得
而我们想要一个最小值，即是该函数值等于0的时候的点，通过求解
得到下列方程：
将其带入代码
# set c to the minimumc = (y - m*x).mean()print(c)
之后同理，可以得到m*
所以有
m = ((y - c)*x).sum()/(x**2).sum()print(m)
然后我们需要创建一个测试点向量，以便将预测结果输入空间来了解拟合程度。
import numpy as npx_test = np.linspace(1890, , 130)[:, None]f_test = m*x_test + c
之后将之前设定的m和c带入，可以看到拟合的结果
plt.plot(x_test, f_test, 'b-')plt.plot(x, y, 'rx')
可以看出拟合的成果并不是很理想。我们尝试给他外面再加一层循环，使其迭代十次。
for i in np.arange(10): m = ((y - c)*x).sum()/(x*x).sum() c = (y-m*x).sum()/y.shape[0]print(m)print(c)f_test = m*x_test + cplt.plot(x_test, f_test, 'b-')plt.plot(x, y, 'rx')最后拟合的效果还是不理想，目标函数的值依旧很大。所以我们需要用目标函数的变化率为基础，直到目标函数找到最小值为止。E1=100000 #尽量大的数E2=((y-m*x-c)**2).sum()while ((E1-E2)/E1)>0.0001: E1=E2 m = ((y-c)*x).sum()/(x*x).sum() c = (y-m*x).sum()/y.shape[0] E2=((y-m*x-c)**2).sum()print(E2)f_test = m * x_test + cplt.plot(x_test, f_test, 'b-')plt.plot(x, y, 'rx')plt.show()
这也就是基于梯度下降算法所优化目标函数而作出的拟合直线。可以看到这是一条符合规律的直线。
二、特解法
基于历史的马拉松数据，设立拟合直线
，所以可设向量矩阵
，
# define the vector ww = np.zeros(shape=(2, 1))w[0] = mw[1] = c
这给了我们原始运算和向量空间中的运算之间的等价性。一般来说，我们现在正转向一个系统，在这个系统中，我们的每一个预测都是由一个内积给出的。当我们想在线性代数中表示一个线性积时，我们倾向于用转置运算来表示，所以我们有
，可以写作
之后我们假设每个数据点xi，都是通过在原始向量上附加1来编写的，从而有
然后我们将所有数据归到矩阵X，有
在代码中可以写成
X = np.hstack((np.ones_like(x), x))print(X)
则有目标函数
其中
，那么
对应的，可以yi被表示成
误差公式原式为，现在可写为，其中
，故代码有
f = np.dot(X, w) # np.dot does matrix multiplication in pythonresid = (y-f)E = np.dot(resid.T, resid) # matrix multiplication on a single vector is equivalent to a dot product.print("Error function is:", E)
为了优化目标函数，需要将目标函数化简，即
有向量矩阵公式
而其中a=c=y以及b=d=Xw，并将其带入原式替换，得到
对w求偏导，可得
为了最小化目标函数，所以会在为0时，有
从而得到
故而拟合的结果如下图所示。可以看出这是一条符合规律的曲线。
w = np.linalg.solve(np.dot(X.T, X), np.dot(X.T, y))print(w)m = w[1]; c=w[0]f_test = m*x_test + cprint(m)print(c)plt.plot(x_test, f_test, 'b-')plt.plot(x, y, 'rx')
三、多元线性回归
以上方法简单的讲述了在一元线性回归中的应用。而此方法也可以建立在多元线性回归上，因为矩阵演算没有指定向量x的长度。首先，我们导入电影数据，并打印电影数据中的列。
import podsdata = pods.datasets.movie_body_count()movies = data['Y']print(', '.join(movies.columns))
现在我们将建立一个基于数字特征的设计矩阵：Year、Body_Count、Length_Minutes，以及IMDB_Rating。我们构建设计矩阵如下：
select_features = ['Year', 'Body_Count', 'Length_Minutes']X = movies.loc[:, select_features]X['Eins'] = 1 # add a column for the offsety = movies[['IMDB_Rating']]
之后根据公式
可得
import pandas as pdimport numpy as npw = pd.DataFrame(data=np.linalg.solve(np.dot(X.T, X), np.dot(X.T, y)), # solve linear regression here index = X.columns, # columns of X become rows of w columns=['regression_coefficient']) # the column of X is the value of regression coefficient
并打印w
w
可以看到IMDB_Rating的值变化相对稳定，与年份和个体数成反比，与电影的时长成正比。而利用矩阵的除法进行w的计算并不是一直稳定的方法，毕竟要把：算出来是要把X的所有元素平方，再除到等式那边去，不仅浪费算力，还会降低数据的精确度。所以便有了QR因式分解。矩阵的QR分解是将其分解为正交矩阵Q，因此，以及一个三角矩阵R。推导公式如下，
代码即结果如下
import scipy as spQ, R = np.linalg.qr(X)w = sp.linalg.solve_triangular(R, np.dot(Q.T, y)) w = pd.DataFrame(w, index=X.columns)w
