最小二乘法擬合
在函數(shù)擬合中丸卷，如果用p表示函數(shù)中需要確定的參數(shù)菌湃，那么目標(biāo)就是找到一組p共缕，使得下面函數(shù)S的值最胁判础：
這種算法稱(chēng)為最小二乘法擬合锯玛。Python的Scipy數(shù)值計(jì)算庫(kù)中的optimize模塊提供了leastsq()函數(shù)，可以對(duì)數(shù)據(jù)進(jìn)行最小二乘擬合計(jì)算。
此處利用該函數(shù)對(duì)一段弧線使用圓方程進(jìn)行了擬合，并通過(guò)Matplotlib模塊進(jìn)行了作圖马僻，程序內(nèi)容如下：
#使用最小二乘法擬合圓方程
import numpy as np
from scipy.optimize import leastsq
from pylab import *
x,y=np.loadtxt('num.dat',delimiter=' ',dtype=float,usecols=(0,1),unpack=True)
#for i in x:
# print i
#for j in y:
# print j
def residuals(p):
a,b,r=p
return r**2-(y-b)**2-(x-a)**2
result=leastsq(residuals,[1,1,1])
a,b,r=result[0]
print "a=",a,"b=",b,"r=",r
plot(x,y,color="red",label="origin",linewidth=2)
yfit=sqrt(r**2-(x-a)**2)+b
plot(x,yfit,"b--",label="fit",linewidth=2)
legend(loc="upper right",frameon=False)
annotate('(x-112.892403261)$^2$+(y+58.8238235027)$^2$=110.123575696$^2$',xy=(40,15))
annotate('R=110.123575696',xy=(100,10))
show()
擬合過(guò)程說(shuō)明
導(dǎo)入模塊
Python的使用中需要導(dǎo)入相應(yīng)的模塊，此處首先用import語(yǔ)句
import numpy as np
from scipy.optimize import leastsq
from pylab import *
分別導(dǎo)入了numpy, leastsq與pylab模塊注服，其中numpy模塊常用用與數(shù)組類(lèi)型的建立韭邓，讀入等過(guò)程。leastsq則為最小二乘法擬合函數(shù)溶弟。pylab是繪圖模塊女淑。
讀入數(shù)據(jù)
接下來(lái)我們需要讀入需要進(jìn)行擬合的數(shù)據(jù)，這里使用了numpy.loadtxt()函數(shù)：
x,y=np.loadtxt('num.dat',delimiter=' ',dtype=float,usecols=(0,1),unpack=True)
其參數(shù)有：
fname：讀入文件名辜御，也可省略參數(shù)直接寫(xiě)'filename'诗力。
dtype：讀入后數(shù)據(jù)的存放類(lèi)型，默認(rèn)為float我抠。
comments：注釋符號(hào)，讀入時(shí)忽略以該符號(hào)開(kāi)頭的行袜茧，默認(rèn)為'#'菜拓。
delimiter：分隔符，存在多列數(shù)據(jù)時(shí)以該符號(hào)進(jìn)行分割笛厦，默認(rèn)為空格纳鼎。
converters：將列號(hào)映射到所需值，可以為遺失的數(shù)據(jù)提供默認(rèn)值裳凸，默認(rèn)為不開(kāi)啟贱鄙。
skiprows：讀入時(shí)跳過(guò)前多少行，默認(rèn)為0姨谷。
usecols：讀取第多少列逗宁，從0開(kāi)始計(jì)數(shù)，可跳躍讀取梦湘，如usecols=(0,3,8)瞎颗。1.11.0版更新，取單獨(dú)一列時(shí)可使用usecols=3捌议，并非必須用元組的形式哼拔。
unpack：布爾值，如果為真瓣颅，則對(duì)返回的數(shù)組進(jìn)行置換倦逐，以便可以使用x,y,z = loadtxt(…)來(lái)解壓縮數(shù)組。當(dāng)與結(jié)構(gòu)化數(shù)據(jù)類(lèi)型一起使用時(shí)宫补，將為每個(gè)字段返回?cái)?shù)組檬姥。默認(rèn)為False曾我。
ndmin：返回的數(shù)組將至少具有ndmin維度。否則穿铆，一維軸將被壓縮您单。合法值:0(默認(rèn)值)、1或2荞雏。1.6.0版本中的新功能虐秦。
encoding：用于解碼。1.14.0版本中的新功能凤优。（不是很理解英文原文的意思）
max_rows：skiprows后讀入的最大行數(shù)悦陋，默認(rèn)讀取全部。1.16.0版本中的新功能筑辨。
進(jìn)行擬合
進(jìn)行擬合時(shí)俺驶，首先我們需要定義一個(gè)目標(biāo)函數(shù)。對(duì)于圓的方程棍辕，我們需要圓心坐標(biāo)(a,b)以及半徑r三個(gè)參數(shù)暮现，方便起見(jiàn)用p來(lái)存儲(chǔ)：
def residuals(p):
a,b,r=p
return r**2-(y-b)**2-(x-a)**2
緊接著就可以進(jìn)行擬合了，leastsq()函數(shù)需要至少提供擬合的函數(shù)名與參數(shù)的初始值：
result=leastsq(residuals,[1,1,1])
a,b,r=result[0]
print "a=",a,"b=",b,"r=",r
返回的結(jié)果為一數(shù)組楚昭，分別為擬合得到的參數(shù)與其誤差值等栖袋，這里只取擬合參數(shù)值。
leastsq()的參數(shù)具體有：
func：函數(shù)名抚太，該函數(shù)應(yīng)至少具有一個(gè)參數(shù)并返回多個(gè)浮點(diǎn)數(shù)塘幅，不能返回NaNs
x0：需擬合參數(shù)的初始值
Dfun：用行間導(dǎo)數(shù)計(jì)算func的雅可比行列的函數(shù)或方法。如果為空尿贫，則估計(jì)雅可比行列式电媳。
full_output：布爾值，非零以返回所有可選輸出庆亡。
col_deriv：布爾值匾乓，非零以指定雅可比函數(shù)計(jì)算列中的導(dǎo)數(shù)（更快，因?yàn)闆](méi)有轉(zhuǎn)置操作）又谋。
ftol：平方和中所需的相對(duì)誤差钝尸。
xtol：近似解中所需的相對(duì)誤差。
gtol：函數(shù)向量和雅可比行列之間所需的正交性搂根。
maxfev：對(duì)函數(shù)的最大調(diào)用次數(shù)珍促。如果提供了Dfun，則默認(rèn)maxfev為100*(N+1)剩愧，其中N是x0中元素的數(shù)量猪叙，否則默認(rèn)maxfev為200*(N+1)。
epsfn：用于確定雅可比行列式的前向誤差合適步長(zhǎng)的變量（對(duì)于Dfun=None）。通常穴翩，實(shí)際步長(zhǎng)將為sqrt(epsfcn)*x如果epsfcn小于機(jī)器精度犬第，則假定相對(duì)誤差為機(jī)器精度的數(shù)量級(jí)。
factor：確定初始步長(zhǎng)的參數(shù)(factor * || diag * x||)芒帕，應(yīng)在(0.1,100)的區(qū)間內(nèi)歉嗓。
diag：N個(gè)正條目，用作變量的比例因子背蟆。
輸出選項(xiàng)有：
x：解（或擬合不成功時(shí)的最后一次迭代的結(jié)果）鉴分。
cov_x：使用fjac和ipvt可選輸出來(lái)構(gòu)建解決方案周?chē)趴杀鹊墓烙?jì)值。如果遇到奇異矩陣則表示無(wú)（表示在某個(gè)方向上的曲率非常平坦）带膀。該矩陣必須乘以殘差方差志珍，得到參數(shù)估計(jì)的協(xié)方差.詳見(jiàn)curve_fit。
infodict：密鑰s的可選輸出字典：
- nfev：函數(shù)調(diào)用的次數(shù)垛叨。
- fvec：在輸出處評(píng)估函數(shù)伦糯。
- fjac：列存儲(chǔ)最終近似雅可比矩陣的QR分解的R矩陣的排列。 與
ipvt一起嗽元，估計(jì)的協(xié)方差的可能近似敛纲。- ipvt：長(zhǎng)度為N的整數(shù)數(shù)組，它定義了一個(gè)置換矩陣p剂癌，使得
fjac*p=q*r淤翔，其中r是上三角形，對(duì)角元素的幅度不增大珍手。 p的列j是單位矩陣的列ipvt(j)。- qtf：向量（轉(zhuǎn)置
(q)*fvec）辞做。
mesg：一個(gè)字符串消息琳要，提供有關(guān)失敗原因的信息。
ies：整數(shù)標(biāo)志秤茅。如果它等于1,2,3或4稚补，則發(fā)現(xiàn)解決方案。否則框喳，找不到解決方案课幕。在任何一種情況下，可選的輸出變量mesg都會(huì)提供更多信息五垮。
對(duì)結(jié)果作圖
最后我們可以將原數(shù)據(jù)與擬合結(jié)果一同做成線狀圖乍惊，可采用pylab.plot()函數(shù)：
plot(x,y,color="red",label="origin",linewidth=2)
plot(x,yfit,"b--",label="fit",linewidth=2)
legend(loc="upper right",frameon=False)
annotate('R=110.123575696',xy=(100,10))
show()
pylab.plot()函數(shù)需提供兩列數(shù)組作為輸入，其他參數(shù)可調(diào)控線條顏色放仗，形狀润绎，粗細(xì)以及對(duì)應(yīng)名稱(chēng)等性質(zhì)。視需求而定，此處不做詳解莉撇。
pylab.legend()函數(shù)可以調(diào)控圖像標(biāo)簽的位置呢蛤，有無(wú)邊框等性質(zhì)。
pylab.annotate()函數(shù)設(shè)置注釋?zhuān)柚辽偬峁┳⑨寖?nèi)容與放置位置坐標(biāo)的參數(shù)棍郎。
pylab.show()函數(shù)用于顯示圖像其障。
最終結(jié)果如下圖所示：
