1. 程式人生 > 其它 >Scipy中最小二乘函式leastsq()

Scipy中最小二乘函式leastsq()

概述

最小二乘法在某種程度上無異於機器學習中基礎中的基礎,且具有相當重要的地位。
optimize模組中提供了很多數值優化演算法,其中,最小二乘法可以說是最經典的數值優化技術了, 通過最小化誤差的平方來尋找最符合資料的曲線。在optimize模組中,使用leastsq()函式可以很快速地使用最小二乘法對資料進行擬合。
比如,有一個未知係數的二元二次函式f(x,y)=w0 x^2 + w1 y^2+w2xy+w3x+w4y+w5,這裡w0~w5為未知的引數,為了確定下來這些引數,將會給定一些樣本點(xi,yi,f(xi,yi)),然後通過調整這些引數,找到這樣一組w0 ~ w5,使得這些所有的樣本點距離函式f(x,y)的距離平方之和最小。
首先看看函式的呼叫格式:

scipy.optimize.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)

引數還是非常多的,一般來說,我們只需要前三個引數就夠了他們的作用分別是:

  • func:誤差函式
  • x0:表示函式的引數
  • args()表示資料點

leastsq,它可以省去中間那些具體的求解步驟,只需要輸入一系列樣本點,給出待求函式的基本形狀(如我剛才所說,二元二次函式就是一種形狀——f(x,y)=w0 x^2 + w1y^2 + w2xy + w3x + w4y + w5,在形狀給定後,我們只需要求解相應的係數w0~w6),即可得到相應的引數。至於中間到底是怎麼求的,這一部分內容就像一個黑箱一樣。(類似與梯度下降法)

例子

函式形為y=kx+b

Xi=np.array([8.19,2.72,6.39,8.71,4.7,2.66,3.78])
Yi=np.array([7.01,2.78,6.47,6.71,4.1,4.23,4.05])

則使用leastsq函式求解其擬合直線的程式碼如下:

###最小二乘法試驗###
import numpy as np
from scipy.optimize import leastsq

###取樣點(Xi,Yi)###
Xi=np.array([8.19,2.72,6.39,8.71,4.7,2.66,3.78])
Yi=np.array([7.01,2.78,6.47,6.71,4.1,4.23,4.05])

###需要擬合的函式func及誤差error###
def func(p,x):
    k,b=p
    return k*x+b

def error(p,x,y,s):
    print s
    return func(p,x)-y #x、y都是列表,故返回值也是個列表

#TEST
p0=[100,2]
#print( error(p0,Xi,Yi) )

###主函式從此開始###
s="Test the number of iteration" #試驗最小二乘法函式leastsq得呼叫幾次error函式才能找到使得均方誤差之和最小的k、b
Para=leastsq(error,p0,args=(Xi,Yi,s)) #把error函式中除了p以外的引數打包到args中
k,b=Para[0]
print"k=",k,'\n',"b=",b

###繪圖,看擬合效果###
import matplotlib.pyplot as plt

plt.figure(figsize=(8,6))
plt.scatter(Xi,Yi,color="red",label="Sample Point",linewidth=3) #畫樣本點
x=np.linspace(0,10,1000)
y=k*x+b
plt.plot(x,y,color="orange",label="Fitting Line",linewidth=2) #畫擬合直線
plt.legend()
plt.show()

1、p0裡放的是k、b的初始值,這個值可以隨意指定。往後隨著迭代次數增加,k、b將會不斷變化,使得error函式的值越來越小。

2、func函式裡指出了待擬合函式的函式形狀。

3、error函式為誤差函式,我們的目標就是不斷調整k和b使得error不斷減小。這裡的error函式和神經網路中常說的cost函式實際上是一回事,只不過這裡更簡單些而已。

4、必須注意一點,傳入leastsq函式的引數可以有多個,但必須把引數的初始值p0和其它引數分開放。其它引數應打包到args中。

5、leastsq的返回值是一個tuple,它裡面有兩個元素,第一個元素是k、b的求解結果,第二個元素我暫時也不知道是什麼意思,先留下來。

其擬合效果圖如下:

函式形為y=ax^2+bx+c

這一次我們給出函式形y=ax^2+bx+c。這種情況下,待確定的引數有3個:a,b和c。

此時給出7個樣本點如下:

Xi=np.array([0,1,2,3,-1,-2,-3])
Yi=np.array([-1.21,1.9,3.2,10.3,2.2,3.71,8.7])
###最小二乘法試驗###
import numpy as np
from scipy.optimize import leastsq

###取樣點(Xi,Yi)###
Xi=np.array([0,1,2,3,-1,-2,-3])
Yi=np.array([-1.21,1.9,3.2,10.3,2.2,3.71,8.7])

###需要擬合的函式func及誤差error###
def func(p,x):
    a,b,c=p
    return a*x**2+b*x+c

def error(p,x,y,s):
    print s
    return func(p,x)-y #x、y都是列表,故返回值也是個列表

#TEST
p0=[5,2,10]
#print( error(p0,Xi,Yi) )

###主函式從此開始###
s="Test the number of iteration" #試驗最小二乘法函式leastsq得呼叫幾次error函式才能找到使得均方誤差之和最小的a~c
Para=leastsq(error,p0,args=(Xi,Yi,s)) #把error函式中除了p以外的引數打包到args中
a,b,c=Para[0]
print"a=",a,'\n',"b=",b,"c=",c

###繪圖,看擬合效果###
import matplotlib.pyplot as plt

plt.figure(figsize=(8,6))
plt.scatter(Xi,Yi,color="red",label="Sample Point",linewidth=3) #畫樣本點
x=np.linspace(-5,5,1000)
y=a*x**2+b*x+c
plt.plot(x,y,color="orange",label="Fitting Curve",linewidth=2) #畫擬合曲線
plt.legend()
plt.show()