【統(tǒng)計(jì)學(xué)習(xí)方法】 從零開始，用python實(shí)現(xiàn)最小二乘法
點(diǎn)擊上方“公眾號”可訂閱哦！
本文主要介紹使用Numpy中的poly1d()多項(xiàng)式函數(shù)和Scipy的optimize模塊的leastsq()函數(shù)，分別使用它們作為擬合函數(shù)和最小二乘法實(shí)現(xiàn)對數(shù)據(jù)的擬合。
1
●
poly1d()函數(shù)
np.poly1d(c_or_r, r=False, variable=None)
該函數(shù)包括三個(gè)參數(shù)，下面分別介紹這三個(gè)參數(shù)。
c_or_r：類似數(shù)組，多項(xiàng)式的系數(shù)，以冪次遞減，或者當(dāng)?shù)诙€(gè)參數(shù)的值為True，即多項(xiàng)式的根（多項(xiàng)式求值為0的值）。
``poly1d（[1，2，3]）``
返回一個(gè)表示以下內(nèi)容的對象：
math：`x ^ 2 + 2x + 3`，
而``poly1d（[1,2,3]，True）``
返回一個(gè)代表：
math：`（x-1）（x-2）（x-3）= x ^ 3-6x ^ 2 + 11x -6`。
r=False：布爾型，可選如果為True，則`c_or_r`指定多項(xiàng)式的根。默認(rèn)為False。
variable：str，可選將打印p時(shí)使用的變量從x更改為variable。
p = np.poly1d([1, 2, 3])print(np.poly1d(p))# output21 x + 2 x + 3?p(0.5)?#?output?計(jì)算當(dāng)x為0.5時(shí)函數(shù)值4.25p.r?#?output?計(jì)算函數(shù)根array([-1.+1.41421356j, -1.-1.41421356j])p.c#?output 現(xiàn)實(shí)系數(shù)array([1, 2, 3])p[1]#?output?顯示多項(xiàng)式中第k次冪的系數(shù)2p * p# output 多項(xiàng)式相乘poly1d([ 1, 4, 10, 12, 9])p = np.poly1d([1,2,3], variable='z')print(p)# output21 z + 2 z + 3# 從其根構(gòu)造一個(gè)多項(xiàng)式：np.poly1d([1, 2], True)# outputpoly1d([ 1, -3, 2])
2
●
leastsq()函數(shù)
首先來看leastsq()函數(shù)的參數(shù)，
leastsq(func,x0,args=(),Dfun=None,full_output=0,col_deriv=0,ftol=1.49012e-08,xtol=1.49012e-08,gtol=0.0,maxfev=0,epsfcn=None,factor=100,diag=None,)
參數(shù)還是非常多的，一般來說，我們只需要前三個(gè)參數(shù)就夠了，
他們的作用分別是：
func：誤差函數(shù)
x0：表示函數(shù)的參數(shù)
args=()：表示數(shù)據(jù)點(diǎn)
3
●
使用最小二乘法進(jìn)行數(shù)據(jù)擬合
import numpy as npimport scipy as spfrom scipy.optimize import leastsqimport matplotlib.pyplot as plt%matplotlib inlinedef func(x):return 2*np.sin(2*np.pi*x)def residuals(p, x, y):fun = np.poly1d(p) # poly1d（）函數(shù)可以按照輸入的列表p返回一個(gè)多項(xiàng)式函數(shù)return y - fun(x) # 返回真實(shí)值 與我們擬合的曲線上對應(yīng)的值的差# 擬合函數(shù)def fitting(p):pars = np.random.rand(p+1) # 生成p+1個(gè)隨機(jī)數(shù)的列表，這樣poly1d函數(shù)返回的多項(xiàng)式次數(shù)就是pr = leastsq(residuals, pars, args=(X, Y)) # 三個(gè)參數(shù)：誤差函數(shù)、函數(shù)參數(shù)列表、數(shù)據(jù)點(diǎn)return r# 要進(jìn)行擬合的數(shù)據(jù)點(diǎn)X = np.linspace(0, 1, 10)Y = [np.random.normal(0, 0.1)+num for num in func(X)] # 添加噪聲# 方便繪制曲線，所以創(chuàng)建多一些點(diǎn)x_ = np.linspace(0, 1, 100)y_ = func(x_)?fit_pars = fitting(3)[0] # 注意返回值中的第一行才是擬合曲線的參數(shù)列表plt.plot(x_, y_, label='real line')plt.scatter(X, Y, label='real points')plt.plot(x_, np.poly1d(fit_pars)(x_), label='fitting line')plt.legend()plt.show()??
?END
掃碼關(guān)注
微信號｜sdxx_rmbj