天天看點

python中matplotlib實作最小二乘法拟合的過程詳解前言一、最小二乘法拟合直線

前言

最小二乘法Least Square Method,做為分類回歸算法的基礎,有着悠久的曆史(由馬裡·勒讓德于1806年提出)。它通過最小化誤差的平方和尋找資料的最佳函數比對。利用最小二乘法可以簡便地求得未知的資料,并使得這些求得的資料與實際資料之間誤差的平方和為最小。最小二乘法還可用于曲線拟合。其他一些優化問題也可通過最小化能量或最大化熵用最小二乘法來表達。

下面這篇文章主要跟大家介紹了關于python中matplotlib實作最小二乘法拟合的相關内容,下面話不多說,來一起看看詳細的介紹:

一、最小二乘法拟合直線

1、拟合直線

生成樣本點

首先,我們在直線 y = 3 + 5x 附近生成服從正态分布的随機點,作為拟合直線的樣本點。

import numpy as np  
import matplotlib.pyplot as plt 
  
# 在直線 y = 3 + 5x 附近生成随機點 
X = np.arange(0, 5, 0.1)  
Z = [3 + 5 * x for x in X]  
Y = [np.random.normal(z, 0.5) for z in Z] 
  
plt.plot(X, Y, 'ro')  
plt.show()
           

樣本點如圖所示:

python中matplotlib實作最小二乘法拟合的過程詳解前言一、最小二乘法拟合直線

拟合直線

設 y = a0 + a1*x,我們利用最小二乘法的正則方程組來求解未知系數 a0 與 a1。

python中matplotlib實作最小二乘法拟合的過程詳解前言一、最小二乘法拟合直線

numpy 的 linalg 子產品中有一個 solve 函數,它可以根據方程組的系數矩陣和方程右端構成的向量來求解未知量。

def linear_regression(x, y):  
 N = len(x) 
 sumx = sum(x) 
 sumy = sum(y) 
 sumx2 = sum(x**2) 
 sumxy = sum(x*y) 
  
 A = np.mat([[N, sumx], [sumx, sumx2]]) 
 b = np.array([sumy, sumxy]) 
  
 return np.linalg.solve(A, b) 
  
a0, a1 = linear_regression(X, Y)
           

繪制直線

此時,我們已經得到了拟合後的直線方程系數 a0 和 a1。接下來,我們繪制出這條直線,并與樣本點做對比。

# 生成拟合直線的繪制點 
_X = [0, 5]  
_Y = [a0 + a1 * x for x in _X] 
  
plt.plot(X, Y, 'ro', _X, _Y, 'b', linewidth=2)  
plt.title("y = {} + {}x".format(a0, a1))  
plt.show()
           

拟合效果如下:

python中matplotlib實作最小二乘法拟合的過程詳解前言一、最小二乘法拟合直線

2、最小二乘法拟合曲線

生成樣本點

與生成直線樣本點相同,我們在曲線 y = 2 + 3x + 4x^2 附近生成服從正态分布的随機點,作為拟合曲線的樣本點。

import numpy as np  
import matplotlib.pyplot as plt 
  
# y = 2 + 3x + 4x^2 
X = np.arange(0, 5, 0.1)  
Z = [2 + 3 * x + 4 * x ** 2 for x in X]  
Y = np.array([np.random.normal(z,3) for z in Z]) 
  
plt.plot(X, Y, 'ro')  
plt.show()
           

樣本點如圖所示:

python中matplotlib實作最小二乘法拟合的過程詳解前言一、最小二乘法拟合直線

拟合曲線

設該曲線的方程為 y = a0 + a1x + a2x^2,同樣,我們通過正則方程組來求解未知量 a0、a1 和 a2。

python中matplotlib實作最小二乘法拟合的過程詳解前言一、最小二乘法拟合直線
#求解多項式數
def quadratic_fitting(X, Y, m):
    N = len(X)
    A = []
    b = []
    # 生成系數矩陣A,b
    for i in range(m+1):
        b.append(sum(X ** i * Y)) # 計算方程組的右端向量b
        # 計算目前方程中的每一個系數
        a = []
        for j in range(m+1):
            a.append(sum(X ** (i + j)))
        A.append(a)
    return np.linalg.solve(A, b) #求解多項式系數
    
 a0, a1, a2 = quadratic_fitting(X, Y, 2)
           

繪制曲線

繪制曲線

我們根據求得的曲線方程,繪制出曲線的圖像。
# 生成拟合曲線的繪制點 
_X = np.arange(0, 5, 0.1)  
_Y = np.array([a0 + a1*x + a2*x**2 for x in _X]) 
  
plt.plot(X, Y, 'ro', _X, _Y, 'b', linewidth=2)  
plt.title("y = {} + {}x + {}$x^2$ ".format(a0, a1, a2))  
plt.show()
           

拟合效果如下:

python中matplotlib實作最小二乘法拟合的過程詳解前言一、最小二乘法拟合直線

我們根據求得的曲線方程,繪制出曲線的圖像。

3、完成代碼

def Lleast_squares_fitting():
    pass


def cubic_polynomial():
    pass #quadratic polynomial


import numpy as np
import matplotlib.pyplot as plt


def linear_fitting(x, y): #linear_regression
    N = len(x)
    sumx = sum(x)
    sumy = sum(y)
    sumx2 = sum(x ** 2)
    sumxy = sum(x * y)

    A = np.mat([[N, sumx], [sumx, sumx2]])
    b = np.array([sumy, sumxy])

    return np.linalg.solve(A, b)

def line_polynomial():
    #首先,我們在直線 y = 3 + 5x 附近生成服從正态分布的随機點,作為拟合直線的樣本點。
    X = np.arange(0, 5, 0.1)
    Z = [3 + 5 * x for x in X]
    Y = [np.random.normal(z, 0.5) for z in Z]

    a0, a1 = linear_fitting(X, Y) ## y = a0 + a1*x 進行拟合
    _X = [0, 5]
    _Y = [a0 + a1 * x for x in _X]

    plt.plot(X, Y, 'ro', _X, _Y, 'b', linewidth=2)
    plt.title("y = {} + {}x".format(a0, a1))

#求解多項式數
def quadratic_fitting(X, Y, m):
    N = len(X)
    A = []
    b = []
    # 生成系數矩陣A,b
    for i in range(m+1):
        b.append(sum(X ** i * Y)) # 計算方程組的右端向量b
        # 計算目前方程中的每一個系數
        a = []
        for j in range(m+1):
            a.append(sum(X ** (i + j)))
        A.append(a)
    return np.linalg.solve(A, b) #求解多項式系數


def quadratic_polynomial():
    #生成樣本點
    # 與生成直線樣本點相同,我們在曲線 y = 2 + 3x + 4x^2 附近生成服從正态分布的随機點,作為拟合曲線的樣本點。
    X = np.arange(0, 5, 0.1)
    Z = [2 + 3 * x + 4 * x ** 2 for x in X]
    Y = np.array([np.random.normal(z, 3) for z in Z])
    # 使用 y = a0 + a1*x + a2*x^2 進行拟合
    a0, a1, a2 = quadratic_fitting(X, Y, 2)

    _X = np.arange(0, 5, 0.1)
    _Y = np.array([a0 + a1 * x + a2 * x ** 2 for x in _X])
    plt.plot(X, Y, 'ro', _X, _Y, 'b', linewidth=2)
    plt.title("y = {} + {}x + {}$x^2$ ".format(a0, a1, a2))


if __name__ == '__main__':
    # line_polynomial()
    quadratic_polynomial()
    plt.show()