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()