劉海靜 劉暢 石慧升 張露予 趙歡
【摘要】本文主要詳細(xì)介紹了如何應(yīng)用Python語言的SciPy庫中的兩種方法實(shí)現(xiàn)一次函數(shù)模型的最小二乘法擬合優(yōu)化計算,通過程序分析及結(jié)果比較得出每種方法的效率及精確性。
【關(guān)鍵詞】SciPy? 最小二乘法? 擬合
1、引言
最小二乘法原理最早出現(xiàn)在Legendre在1805年發(fā)表的《計算彗星軌道的新方法》的附錄中。Legendre認(rèn)為:“賦予誤差的平方和為極小時,則意味著在這些誤差之間建立了一種均衡性,它阻止了極端情形所施加的過度影響。這非常好地適用于揭示最接近真實(shí)情形的系統(tǒng)狀態(tài)”。
2、理論基礎(chǔ)
2.1最小二乘法理論基礎(chǔ)
最小二乘法主要應(yīng)用于解決連續(xù)或離散函數(shù)模型的最優(yōu)解問題,是測量工作以及其它科學(xué)工程領(lǐng)域中,應(yīng)用最早也是最廣泛的算法之一。最小二乘法是一種最常見的曲線擬合數(shù)學(xué)優(yōu)化方法,其原理是將復(fù)雜函數(shù)變?yōu)楹唵魏瘮?shù)或者將離散函數(shù)轉(zhuǎn)變?yōu)檫B續(xù)函數(shù),其意義在于,簡化函數(shù)或數(shù)學(xué)模型,將原先不易求解和解析困難的問題轉(zhuǎn)化為容易求解和容易析解的問題。其通常的解析過程需要首先在一個連續(xù)或者離散的原函數(shù)上采樣若干樣本點(diǎn),接著用若干個簡單的函數(shù)與采樣樣本點(diǎn)進(jìn)行擬合,并求擬合函數(shù)與采樣點(diǎn)之間的方差,循環(huán)重復(fù)這一過程,直到樣本中出現(xiàn)一個最方差,這條擬合曲線就是通過最小二乘法得到的最優(yōu)擬合曲線。
在生產(chǎn)和實(shí)踐過程中,經(jīng)常會遇到利用一組觀測數(shù)據(jù)通過最小二乘法來估計某些未知參數(shù)的問題。最小二乘法在計算過程中存在一定的誤差,并且使用人工計算的方法計算量非常大,但隨著近代計算機(jī)技術(shù)的不斷發(fā)展,得益于計算機(jī)強(qiáng)大的計算能力,最小二乘法計算量大的問題得到了很好的解決,因此最小二乘法已經(jīng)成為了工程和應(yīng)用數(shù)學(xué)領(lǐng)域最常用的函數(shù)擬合優(yōu)化方法。
2.2 SciPy基礎(chǔ)
SciPy起源于2001年,是一個面向Python語言的開源科學(xué)計算數(shù)據(jù)包,常用用于數(shù)學(xué)、科學(xué)、工程領(lǐng)域,可以處理插值、積分、優(yōu)化、整合以及線性代數(shù)模塊、傅里葉變換、圖像處理、常微分方程數(shù)值解的求解、信號處理等各種數(shù)值計算的問題。SciPy庫是構(gòu)建在NumPy庫的基礎(chǔ)之上的,SciPy提供了許多操作NumPy數(shù)組的函數(shù),可以有效得計算NumPy庫中的矩陣和數(shù)組,使得NumPy庫和Scipy庫協(xié)同工作,高效得解決數(shù)據(jù)計算和數(shù)據(jù)處理問題。
scipy.sparse模塊構(gòu)建了7種稀疏矩陣數(shù)據(jù)結(jié)構(gòu),7種結(jié)構(gòu)在計算機(jī)中的以稀疏格式儲存,使得科學(xué)計算和數(shù)據(jù)處理有效得節(jié)省了空間,提高了運(yùn)行效率。其中最常用的CSR與CSC這兩種方法提供了快速的主軸索引方式和稀疏模式下的快速矩陣向量乘法,大大提高了矩陣的運(yùn)算速度,從而在依托于SciPy庫的科學(xué)計算中得到了廣泛的應(yīng)用。scipy.linalg模塊提供了標(biāo)準(zhǔn)的線性代數(shù)運(yùn)算,同時還加入了NumPy中沒有的高級函數(shù),使用BLAS/LAPACK庫編譯。scipy.optimize模塊提供了非線性方程求解,函數(shù)最小值,曲線擬合和路徑尋找等數(shù)學(xué)優(yōu)化方法。
3、最小二乘法求解設(shè)計實(shí)現(xiàn)
3.1 Scipy.linalg模擬干擾樣本對一次函數(shù)求解
本文首先采用scipy.linalg模塊。
3.1.1設(shè)計思路。利用numpy.linspace在(-3,3)上建立150個初始數(shù)據(jù),建立理想模型y=5x+3,通過numpy.random模塊中的隨機(jī)函數(shù)模擬μ為0,σ為0.3的正態(tài)分布外界干擾,并添加種子池,保證每次程序運(yùn)行干擾值固定,對初始測試樣本進(jìn)行干擾,通過干擾后的測試樣本構(gòu)建數(shù)學(xué)模型方程組,抽取方程組中的系數(shù)矩陣,并進(jìn)行numpy縱向堆疊,傳入scipy.linalg模塊lstsq函數(shù),由于lstsq函數(shù)傳入?yún)?shù)需要為數(shù)組形式,所以對堆疊之后的結(jié)果進(jìn)行轉(zhuǎn)置,再將轉(zhuǎn)置后的結(jié)果傳入lstsq函數(shù)。
lstsq函數(shù)的返回值共有4個,所求的結(jié)果系數(shù)在第一個返回值中,而第一個返回值為數(shù)組類型,根據(jù)數(shù)學(xué)模型,可以對應(yīng)相應(yīng)的索引號打印顯示出最終的最優(yōu)模型系數(shù),最后通過Matplotlib庫進(jìn)行數(shù)據(jù)可視化輸出干擾后的測試樣本和擬合曲線,并在lstsq函數(shù)求解的首尾加入計時函數(shù),具體流程圖如圖1所示。
3.1.2實(shí)現(xiàn)程序
import numpy as np
from scipy import linalg
import matplotlib.pyplot as plt
import time
plt.rcParams['font.sans-serif'] = ['SimHei']? ? ? ?# 用來正常顯示中文標(biāo)簽
plt.rcParams['axes.unicode_minus'] = False? ? ? # 用來正常顯示負(fù)號
start = time.time()
np.random.seed(0)
x=np.linspace(-3,3,150)? ? ? ? ? ? ? #x的定義域,150個模擬的未擾動測試樣本
y=5*x+3? ? ? ? ? ? ? ? ? ? ? ? ? #建立的理想數(shù)學(xué)模型,預(yù)期值
x1=x+np.random.normal(0,0.3,150)? ? #現(xiàn)實(shí)外界條件干擾,通過random函數(shù)模塊模擬外
界干擾