一、最小二乘拟合直线
生成样本点
首先，在直线 y = 3 + 5*x 附近生成服从正态分布的随机样本点，作为拟合直线的样本点，即实际使用中的观测点数据
1 #最小二乘拟合直线 2 import numpy as np 3 import matplotlib.pyplot as plt
1 #随机点 2 x = np.arange(0,5,0.1) 3 z = [3+5*x for x in x] 4 y = [np.random.normal(z,1) for z in z] 5 #设置显示中文 6 plt.rcParams['font.family'] = 'SimHei' 7 #绘图标题 8 plt.title("这里时标题",fontsize = 20) 9 plt.xlabel("这里是x轴",fontsize=14) 10 plt.ylabel("这里是y轴",fontsize=14) 11 12 plt.plot(x,y,'bp',label="point") 13 plt.plot(x,z,'r',label="line") 14 #这是图例 15 plt.legend(prop = {"family":"Times New Roman","size":14}) 16 #plt.savefig() #输出图像 17 plt.show()
如图所示：
拟合直线
设 y = a0 + a1*x ，利用最小二乘正则方程组求解未知系数 a0与a1。
numpy 的 linalg 模块中有一个 solve 函数，它可以根据方程组的系数矩阵和方程右端构成的向量来求解未知量。
1 def linear_regression(x,y): 2 N = len(x) 3 sumx = sum(x) 4 sumy = sum(y) 5 sumx2 = sum(x**2) 6 sumxy = sum(x*y) 7 A = np.mat([[N,sumx],[sumx,sumx2]]) 8 b = np.array([sumy,sumxy]) 9 10 return np.linalg.solve(A,b) 11 a0,a1 = linear_regression(x,y)
此时，我们已经得到了拟合后的直线方程系数 a0 和 a1。接下来，我们绘制出这条直线，并与样本点做对比
1 _x = [0,5] 2 _y = [a0 + a1*x for x in _x] 3 4 #设置显示中文 5 plt.rcParams['font.family'] = 'SimHei' 6 #绘图标题 7 plt.xlabel("这里是x轴",fontsize=14) 8 plt.ylabel("这里是y轴",fontsize=14) 9 plt.plot(x,y,'bo',label='point') 10 plt.plot(_x,_y,'r',label="line") 11 plt.title("y = {} + {}x".format(a0,a1)) 12 #这是图例 13 plt.legend(prop = {"family":"Times New Roman","size":14}) 14 plt.show()
二、最小二乘拟合曲线
与生成直线样本点相同，我们在曲线 y = 2 + 3x + 4x^2 附近生成服从正态分布的随机点，作为拟合曲线的样本点。
1 import numpy as np 2 import matplotlib.pyplot as plt 3 # y = 2+3x+4x^2 4 x = np.arange(0,5,0.1) 5 z = [2+3*x+4*x**2 for x in x] 6 y = np.array([np.random.normal(z,3) for z in z]) 7 8 plt.plot(x,y,'ro') 9 plt.show()
拟合曲线
设该曲线的方程为 y = a0 + a1*x + a2*x^2，同样，我们通过正则方程组来求解未知量 a0、a1 和 a2。
1 # 生成系数矩阵A 2 def gen_coefficient_matrix(X): 3 N = len(X) 4 m = 3 5 A = [] 6 # 计算每一个方程的系数 7 for i in range(m): 8 a = [] 9 # 计算当前方程中的每一个系数 10 for j in range(m): 11 a.append(sum(X ** (i+j))) 12 A.append(a) 13 return A 14 15 # 计算方程组的右端向量b 16 def gen_right_vector(X, Y): 17 N=len(X) 18 m=3 19 b=[] 20 for i in range(m): 21 b.append(sum(X**i*y)) 22 return b 23 24 A = gen_coefficient_matrix(x) 25 b = gen_right_vector(x,y) 26 27 a0, a1, a2 = np.linalg.solve(A, b)
1 _X = np.arange(0, 5, 0.1) 2 _Y = np.array([a0 + a1*x + a2*x**2 for x in _X]) 3 4 plt.plot(x, y, 'ro', _X, _Y, 'b', linewidth=2) 5 plt.title("y = {} + {}x + {}$x^2$ ".format(a0, a1, a2)) 6 plt.show()
【参考】：https://www.jb51.net/article/118406.htm