最小二乘擬合

來自：某小皮

最優化函數庫Optimization

優化是找到最小值或等式的數值解的問題。scipy.optimization子模塊提供函數最小值，曲線擬合和尋找等式的根的有用算法。

 

最小二乘擬合

假設有一組實驗數據（xi, yi），事先知道它們之間應該滿足某函數關系yi = f(xi)，通過這些已知的信息，需要確定函數f的一些參數。例如，如果函數f是線性函數f(x) = kx + b，那么參數k和b就是需要確定的值。

如果用p表示函數中需要確定的參數，那么目標就是找到一組p，使下面的函數s的值最小：

這種算法被稱作最小二乘擬合(Least-square fitting)。

使用leastsq()進行最小二乘擬合計算。leastsq()只需要將計算誤差的函數和待確定參數的初始值傳遞給它即可。

leastsq()函數傳入誤差計算函數和初始值，該初始值將作為誤差計算函數的第一個參數傳入；

計算的結果r是一個包含兩個元素的元組，第一個元素是一個數組，表示擬合后的參數k，b；第二個元素如果等於1，2，3，4中的其中一個整數，則擬合成功，否則將返回mesg。

 

leastsq函數：

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=0.0, factor=100, diag=None, warning=True) 

一般我們只要指定前三個參數：

• func 是我們自己定義的一個計算誤差的函數，
• x0 是計算的初始參數值
• args 是指定func的其他參數

注意：傳入leastsq函數的參數可以有多個，但必須把參數的初始值p0和其它參數分開放。其它參數應打包到args中。

 

最小二乘擬合示例：

使用最小二乘對帶噪聲的正弦波數據進行擬合：

擬合得到的參數雖然和實際的參數有可能完全不同，但是由於正弦函數具有周期性，實際上擬合的結果和實際的函數是一致的。

import numpy as np
from scipy.optimize import leastsq

#待擬合的函數，x是變量，p是參數
def fun(x, p):
    a, b = p
    return a*x + b

#計算真實數據和擬合數據之間的誤差，p是待擬合的參數，x和y分別是對應的真實數據
def residuals(p, x, y):
    return fun(x, p) - y

#一組真實數據，在a=2, b=1的情況下得出
x1 = np.array([1, 2, 3, 4, 5, 6], dtype=float)
y1 = np.array([3, 5, 7, 9, 11, 13], dtype=float)

#調用擬合函數，第一個參數是需要擬合的差值函數，第二個是擬合初始值，第三個是傳入差值函數的其他參數
r = leastsq(residuals, [1, 1], args=(x1, y1))

#打印結果，r[0]存儲的是擬合的結果，r[1]、r[2]代表其他信息
print(r[0])

運行結果：[ 2. 1.]

 

 

import numpy as np
import pylab as pl
from scipy.optimize import leastsq

def func(x, p):
    # 數據擬合所用的函數：A*sin(2*pi*k*x + theta)
    A, k, theta = p
    return A*np.sin(2*np.pi*k*x + theta)
def residuals(p, y, x):
    # 真實數據x, y和擬合數據之間的誤差，p為待擬合參數
    return y - func(x, p)

x = np.linspace(-2*np.pi, 0, 100)

A, k, theta = 10, 0.34, np.pi/6     　　 # 真實數據的函數參數
y0 = func(x, [A, k, theta])          　　# 真實數據
y1 = y0 + 2 * np.random.randn(len(x))　　# 加入噪聲之后的實驗數據

p0 = [7, 0.2, 0]    # 第一次猜測的函數擬合參數

# 調用leastsq進行數據擬合，
# residuals為我們自己定義的需要擬合的差值函數
# p0為擬合初始值，
# args為傳入差值函數的其他參數
plsq = leastsq(residuals, p0, args = (y1, x))

"""
除了初始值之外，還調用了args參數，用於指定residuals中使用到的其他參數（直線擬合是直接使用X,Y的全局變量），
同樣也返回一個元組，第一個元素為擬合后的參數數組；
這里將（y1, x）傳遞給args參數。Leastsq()會將這兩個額外的參數傳遞給residuals()。
因此residuals()有三個參數，p是正弦函數的參數，y和x是表示真實數據的數組。
"""

print(u"真實參數：", [A, k, theta])
print(u"擬合參數", plsq[0])    # 實驗數據擬合后的參數
pl.plot(x, y0, label = u"真實數據")
pl.plot(x, y1, label = u"帶噪聲的數據")
pl.plot(x, func(x, plsq[0]), label = u"擬合數據")
pl.legend()
pl.show()

 

 運行結果：