回归方程在python中的复现#花树P69

正规方程

本周实现花树P69页正规方程在python中的复现。
首先，回归方程的建模公式如下图：

其中y^ 为向量，W为向量，X为矩阵。其中W权值是我们最后所求的，因此我们引入均方误差MSE作为衡量模型最终效果的指标，公式如下：

其中，y^为测试集预测值向量，y为测试机真实值向量。然后对MSE求导，求其倒数为0的点时，W的值：


最后求得式5.12，为正规方程式。用正规方程求得的W，为MSE最小时的W，是最逼近样本点的回归方程。其中比较难的一步是5.10到5.11的求导，这里放上公式：

简单的理解是W.TX.TWX 为一个二次型，而且W.TX.TWX 最后结果是一个标量，因此对W求导，则会求得2X.TX.W。

python中复现回归方程

这里为了方便展示，x我们只取一维x1，因此W也只有一个W1。

# -*- coding: utf-8 -*-
"""
Created on Fri Oct 14 17:07:03 2022

@author: dell
"""

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif']=['SimHei'] #　显示中文
plt.rcParams['axes.unicode_minus']=False #显示正负号


m=10
n=1
w = np.array([[0.5]])
X = np.random.randn(m,n)
#y = np.dot(w.T , X)
y = np.random.randn(m ,n)

# 正规方程
W = np.dot( np.linalg.inv( np.dot( X.T , X)) , np.dot(X.T , y))

# 求解 MSE 和 权值
y_pred = W * X                                                      # 根据正规方程求得的权值
mse_best=((1/m) * np.sqrt( np.linalg.norm(y_pred - y) ** 2))        # 最优权值算出最优mse
mse = []

for WW in np.arange(W-0.3,W + 0.7 ,0.01):                                      # 从0-1筛选权值
    print(WW)
    y_pred = WW * X                                                 # 由0-1的权值得到的y预测值
    mse.append((1/m) * np.sqrt( np.linalg.norm(y_pred - y) ** 2))   #  由0-1的权值得到的mse


plt.figure()
plt.subplot(121)
plt.scatter(X,y.reshape(m))
plt.plot(X,W*X)
plt.title('线性回归示例')
plt.ylabel('y')
plt.xlabel('x1')
plt.subplot(122)
plt.plot(np.arange(W-0.3,W + 0.7 ,0.01) , mse)
plt.scatter(W , mse_best)
plt.title('w的最优解')
plt.xlabel('w1')
plt.ylabel('MSE')
plt.show()

结果图展示：

可以看到，回归方程得到的W是MSE最小的W。