一、最小二乘擬合直線
生成樣本點
首先，在直線 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