线性回归 ：标准方程法的求解 （包含numpy对矩阵加工的代码）

线性回归

生活满意度的回归模型：$$life-satisfaction={\theta }_0+{\theta }_1\times GDP-per-capita$$
这个模型就是输入特征GDP_per_capita的线性函数，${\theta }_0$ 和 ${\theta }_1$是模型的参数。

更为概括地说，线性模型就是对输入特征加权求和，再加上一个我们称为偏置项（也称为截距项）的常数，以此进行预测。
线性回归模型预测公式如下：
$$\hat{y}={\theta }_0+{\theta }_1{x}_1+{\theta }_2{x}_2+\cdots +{\theta }_n{x}_n$$
在此等式中：

• $\hat{y}$是预测值。
• $n$是特征数量。
• $x_i$是第i个特征值。
• ${\theta }_j$是第j个模型参数（从0开始，${\theta }_0$也算）

也可以使用向量化的形式更简洁地表示
$$\hat{y}=h_{\theta }(x)=\theta \times x$$

在此等式中：

• $\theta$是模型的参数向量，其中包含偏差项${\theta }_0$和特征权重${\theta }_1$至
  n。
• $x$是实例的特征向量，包含从$x_0$到$x_n$，$x_0$始终等于1。
• $\theta \times x$是向量$\theta$和$x$的点积，它当然等于${\theta }_0{x}_0+{\theta }_1{x}_1+{\theta }_2{x}_2+\cdots +{\theta }_n{x}_n$
• $h_{\theta }$是假设函数，使用模型参数$\theta$。

在机器学习中，向量通常表示为列向量，是有单一列的二维数组。如果$\theta$和$x$为列向量，则预测为 ，其中${\theta }^T$为$\theta$（行向量而不是列向量）的转置，且${\theta }^T$ $x$为${\theta }^T$和$x$的矩阵乘积。当然这是相同的预测，除了现在是以单一矩阵表示而不是一个标量值。

这就是线性回归模型，我们该怎样训练线性回归模型呢？回想一下，训练模型就是设置模型参数直到模型最拟合训练集的过程。为此，我们首先需要知道怎么测量模型对训练数据的拟合程度是好还是差。我们了解到回归模型最常见的性能指标是均方根误差（RMSE）。

因此，在训练线性回归模型时，你需要找到最小化RMSE的$\theta$值。在实践中，将均方误差（MSE）最小化比最小化RMSE更为简单，二者效果相同（因为使函数最小化的值，同样也使其平方根最小）
MSE公式如下
$$MSE=(X,h_{\theta })=\frac{1}{m}\sum _{i=1}^m({\theta }^T{x}^{(i)}-y^{(i)}{)}^2$$
$y^{(i)}$表示第i个真实值。

标准方程

为了得到使成本函数最小的θ值，有一个闭式解方法——也就是一个直接得出结果的数学方程，即标准方程。
公式：标准方程$$\hat{\theta }=(X^{T}X{)}^{-1}{X}^{T}y$$
这个方程中：

• $\hat{\theta }$是使成本函数最小的θ值。

• $y$是包含$y^{(i)}$到$y^{(m)}$的目标值向量。

import numpy as np
import matplotlib.pyplot as plt
X = 2*np.random.rand(100,1) #生成一个形状为(100 ,1)的数组，从[0,2)中随机取值
y = 4+3*X +np.random.randn(100,1) 

np.random.rand(a,b) 生成一个在[0,1)上随机取值，形状为(a,b)的二维均匀分布数组
np.random.randn(a,b) 生成一个均值为0，标准差为1 ，形状为(a,b)的二维正态分布数组

画图查看分布情况

plt.plot(X,y,"r.")
plt.xlabel("$x_1$",fontsize=18)
plt.ylabel("$y$",fontsize = 18)
plt.axis([0,2,0,13])
plt.savefig("generated_data_plot")
plt.show()

结果如下

现在我们使用标准方程来计算$\hat{\theta }$。使用Numpy 的线性代数模块$(np.linalg)$中的$inv()$函数来对矩阵求逆，并使用$dot()$方法计算矩阵内积

X_b = np.c_[np.ones((100,1)),X] 
theta_best = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y)

$$\hat{\theta }=(X^{T}X{)}^{-1}{X}^{T}y$$
我们实际用来生成数据的函数是$y=4+3x_1+$高斯噪声。
我们要求的是${\theta }_0$ 和${\theta }_1$，${\theta }_0$的系数是1，故需要添加一个1的矩阵。
看看结果如何

array([[4.21509616],
       [2.77011339]])

我们期待的是${\theta }_0=4$和${\theta }_1=3$得到的是${\theta }_0=4.215$，${\theta }_1=2.770$。非常
接近，噪声的存在使其不可能完全还原为原本的函数。

现在可以用$\hat{\theta }$作出预测。

X_new = np.array([[0],[2]])
X_new_b = np.c_[np.ones((2,1)),X_new]
y_predict = X_new_b.dot(theta_best)
y_predict

结果如下：

array([[4.21509616]
 [9.75532293]])

绘制模型的预测结果：

plt.plot(X_new,y_predict,"r-")
plt.plot(X,y,"b.")
plt.axis([0,2,0,15])
plt.show()

结果如下：

使用Scikit-Learn执行线性回归也很简单

from sklearn.linear_model import LinearRegression
lin_reg = LinearRegression()
lin_reg.fit(X,y)
lin_reg.intercept_ #截距 
lin_reg.coef_ #系数
#(array([4.21509616]), array([[2.77011339]]))

LinearRegression类基于scipy.linalg.lstsq（）函数（名称代表“最小二乘”），你可以直接调用它：

theta_best_svd ,residuals,rank, s = np.linalg.lstsq(X_b,y,rcond = 1e-6)
theta_best_svd

函数的返回值是一个元组，包含以下四个元素：

• theta_best_svd：这是最小二乘解，即模型参数的最优估计。在线性回归的上下文中，这些参数对应于自变量（包括截距项，如果添加了的话）的系数。
• residuals：这是残差数组，即实际观测值与模型预测值之间的差异。
• rank：这是设计矩阵X_b的数值秩。它反映了矩阵中独立（非零）的行（或列）的数量。
• s：这是奇异值的数组。这些值是从设计矩阵的奇异值分解中得到的，可以用来评估数值稳定性和模型的复杂度。
  结果如下：

array([[4.21509616],
       [2.77011339]])

此函数计算$\hat{\theta }=X^{+}y$，其中$X^{+}$是$X$的伪逆（具体来说是MoorePenrose逆）。你可以使用$np.linalg.pinv（）$来直接计算这个伪逆：

np.linalg.pinv(X_b).dot(y)

伪逆：$$X^{+}=(X^{T}X{)}^{-1}{X}^T$$