線性模型（二）之多項(xiàng)式擬合
讓客戶滿意是我們工作的目標(biāo)，不斷超越客戶的期望值來(lái)自于我們對(duì)這個(gè)行業(yè)的熱愛。我們立志把好的技術(shù)通過(guò)有效、簡(jiǎn)單的方式提供給客戶，將通過(guò)不懈努力成為客戶在信息化領(lǐng)域值得信任、有價(jià)值的長(zhǎng)期合作伙伴，公司提供的服務(wù)項(xiàng)目有：空間域名、網(wǎng)站空間、營(yíng)銷軟件、網(wǎng)站建設(shè)、察布查爾錫伯網(wǎng)站維護(hù)、網(wǎng)站推廣。
1. 多項(xiàng)式擬合問題
??多項(xiàng)式擬合（polynominal curve fitting）是一種線性模型，模型和擬合參數(shù)的關(guān)系是線性的。多項(xiàng)式擬合的輸入是一維的，即x=xx=x，這是多項(xiàng)式擬合和線性回歸問題的主要區(qū)別之一。
??多項(xiàng)式擬合的目標(biāo)是構(gòu)造輸入xx的MM階多項(xiàng)式函數(shù)，使得該多項(xiàng)式能夠近似表示輸入xx和輸出yy的關(guān)系，雖然實(shí)際上xx和yy的關(guān)系并不一定是多項(xiàng)式，但使用足夠多的階數(shù)，總是可以逼近表示輸入xx和輸出yy的關(guān)系的。
??多項(xiàng)式擬合問題的輸入可以表示如下：
D={(x1,y1),(x2,y2),...,(xi,yi),...,(xN,yN)}xi∈Ryi∈R
D={(x1,y1),(x2,y2),...,(xi,yi),...,(xN,yN)}xi∈Ryi∈R
??目標(biāo)輸出是得到一個(gè)多項(xiàng)式函數(shù)：
f(x)=w1x1+w2x2+wixi+...+wMxM+b=(∑i=1Mwixi)+b
f(x)=w1x1+w2x2+wixi+...+wMxM+b=(∑i=1Mwixi)+b
其中MM表示最高階數(shù)為MM。
??可見在線性擬合的模型中，共包括了(M+1)(M+1)個(gè)參數(shù)，而該模型雖然不是輸入xx的線性函數(shù)，但卻是(M+1)(M+1)個(gè)擬合參數(shù)的線性函數(shù)，所以稱多項(xiàng)式擬合為線性模型。對(duì)于多項(xiàng)式擬合問題，其實(shí)就是要確定這(M+1)(M+1)個(gè)參數(shù)，這里先假設(shè)階數(shù)MM是固定的（MM是一個(gè)超參數(shù)，可以用驗(yàn)證集來(lái)確定MM最優(yōu)的值，詳細(xì)的關(guān)于MM值確定的問題，后面再討論），重點(diǎn)就在于如何求出這(M+1)(M+1)個(gè)參數(shù)的值。
2.優(yōu)化目標(biāo)
??多項(xiàng)式擬合是利用多項(xiàng)式函數(shù)逼近輸入xx和輸出yy的函數(shù)關(guān)系，通過(guò)什么指標(biāo)來(lái)衡量某個(gè)多項(xiàng)式函數(shù)的逼近程度呢？（其實(shí)這就是誤差/損失函數(shù)）。擬合/回歸問題常用的評(píng)價(jià)指標(biāo)是均方誤差（在機(jī)器學(xué)習(xí)中的模型評(píng)估與度量博客中，我進(jìn)行了介紹）。多項(xiàng)式擬合問題也同樣采用該評(píng)價(jià)指標(biāo)，以均方誤差作為誤差/損失函數(shù)，誤差函數(shù)越小，模型越好。
E(w,b)=1N∑i=1N[f(xi)?yi]2
E(w,b)=1N∑i=1N[f(xi)?yi]2
??系數(shù)1N1N是一常數(shù)，對(duì)優(yōu)化結(jié)果無(wú)影響，可以去除，即將均方誤差替換為平方誤差：
E(w,b)=∑i=1N[f(xi)?yi]2
E(w,b)=∑i=1N[f(xi)?yi]2
?? 到這里，就成功把多項(xiàng)式擬合問題變成了最優(yōu)化問題，優(yōu)化問題可表示為：
argminw,bE(w,b)
arg?minw,b?E(w,b)
即需要求得參數(shù){w1,...,wM,b}{w1,...,wM,b}的值，使得E(w,b)E(w,b)最小化。那么如何對(duì)該最優(yōu)化問題求解呢？
3. 優(yōu)化問題求解
3.1 求偏導(dǎo)，聯(lián)立方程求解
?? 直觀的想法是，直接對(duì)所有參數(shù)求偏導(dǎo)，令偏導(dǎo)為0，再聯(lián)立這M+1M+1個(gè)方程求解（因?yàn)楣灿蠱+1M+1個(gè)參數(shù)，故求偏導(dǎo)后也是得到M+1M+1個(gè)方程）。
E(w,b)=∑i=1N[f(xi)?yi]2=∑i=1N[(w1x1i+w2x2i+wixji+...+wMxMi+b)?yi]2
E(w,b)=∑i=1N[f(xi)?yi]2=∑i=1N[(w1xi1+w2xi2+wixij+...+wMxiM+b)?yi]2
利用E(w,b)E(w,b)對(duì)各個(gè)參數(shù)求偏導(dǎo)，如下：
?E(w,b)?wj?E(w,b)?b=2∑i=1N[(w1x1i+w2x2i+wixji+...+wMxMi+b)?yi]xji=2∑i=1N[(w1x1i+w2x2i+wixji+...+wMxMi+b)?yi]
?E(w,b)?wj=2∑i=1N[(w1xi1+w2xi2+wixij+...+wMxiM+b)?yi]xij?E(w,b)?b=2∑i=1N[(w1xi1+w2xi2+wixij+...+wMxiM+b)?yi]
求導(dǎo)之后，將各個(gè)點(diǎn)(xi,yi)(xi,yi)的值帶入偏導(dǎo)公式，聯(lián)立方程求解即可。
??針對(duì)該解法，可以舉個(gè)例子詳細(xì)說(shuō)明，比如有兩個(gè)點(diǎn)(2,3),(5,8)(2,3),(5,8),需要利用二階多項(xiàng)式f(x)=w1x+w2x2+bf(x)=w1x+w2x2+b擬合。求解過(guò)程如下：
該二階多項(xiàng)式對(duì)參數(shù)求偏導(dǎo)得到
?E(w,b)?wj?E(w,b)?b=2∑i=12[(w1x1i+w2x2i+b)?yi]xji=[(w1x1+w2x21+b)?y1]xj1+[(w1x2+w2x22+b)?y2]xj2=2∑i=12[(w1x1i+w2x2i+b)?yi]=[(w1x1+w2x21+b)?y1]+[(w1x2+w2x22+b)?y2]
?E(w,b)?wj=2∑i=12[(w1xi1+w2xi2+b)?yi]xij=[(w1x1+w2x12+b)?y1]x1j+[(w1x2+w2x22+b)?y2]x2j?E(w,b)?b=2∑i=12[(w1xi1+w2xi2+b)?yi]=[(w1x1+w2x12+b)?y1]+[(w1x2+w2x22+b)?y2]
將點(diǎn)(2,3),(5,8)(2,3),(5,8)帶入方程，可以得到3個(gè)方程，
2b+7w1+29w2=117b+29w1+133w2=4629b+133w1+641w2=212
2b+7w1+29w2=117b+29w1+133w2=4629b+133w1+641w2=212
聯(lián)立這三個(gè)方程求解，發(fā)現(xiàn)有無(wú)窮多的解，只能得到3w1+21w2=53w1+21w2=5，這三個(gè)方程是線性相關(guān)的，故沒有唯一解。
??該方法通過(guò)求偏導(dǎo)，再聯(lián)立方程求解，比較復(fù)雜，看著也很不美觀。那么有沒有更加方便的方法呢？
3.2 最小二乘法
?? 其實(shí)求解該最優(yōu)化問題（平方和的最小值）一般會(huì)采用最小二乘法（其實(shí)最小二乘法和求偏導(dǎo)再聯(lián)立方程求解的方法無(wú)本質(zhì)區(qū)別，求偏導(dǎo)也是最小二乘法，只是這里介紹最小二乘的矩陣形式而已）。最小二乘法（least squares），從英文名非常容易想到，該方法就是求解平方和的最小值的方法。
??可以將誤差函數(shù)以矩陣的表示(NN個(gè)點(diǎn)，最高M(jìn)M階)為：
∥Xw?y∥2
‖Xw?y‖2
其中，把偏置bb融合到了參數(shù)ww中，
w={b,w1,w2,...,wM}
w={b,w1,w2,...,wM}
XX則表示輸入矩陣，
??????11...1x1x2...xNx21x22...x2N............xM1xM2...xMN??????
[1x1x12...x1M1x2x22...x2M...............1xNxN2...xNM]
yy則表示標(biāo)注向量，
y={y1,y2,...,yN}T
y={y1,y2,...,yN}T
因此，最優(yōu)化問題可以重新表示為
minw∥Xw?y∥2
minw‖Xw?y‖2
對(duì)其求導(dǎo)，
?∥Xw?y∥2?w=?(Xw?y)T(Xw?y)?w=?(wTXT?yT)(Xw?y)?w=?(wTXTXw?yTXw?wTXTy+yTy)?w
?‖Xw?y‖2?w=?(Xw?y)T(Xw?y)?w=?(wTXT?yT)(Xw?y)?w=?(wTXTXw?yTXw?wTXTy+yTy)?w
在繼續(xù)對(duì)其求導(dǎo)之前，需要先補(bǔ)充一些矩陣求導(dǎo)的先驗(yàn)知識(shí)（常見的一些矩陣求導(dǎo)公式可以參見轉(zhuǎn)載的博客），如下：
?xTa?x=a?ax?x=aT?xTA?x=Ax+ATx
?xTa?x=a?ax?x=aT?xTA?x=Ax+ATx
根據(jù)上面的矩陣求導(dǎo)規(guī)則，繼續(xù)進(jìn)行損失函數(shù)的求導(dǎo)
?∥Xw?y∥2?w=?(wTXTXw?yTXw?wTXTy+yTy)?w=XTXw+(XTX)Tw?(yTX)T?XTy=2XTXw?2XTy
?‖Xw?y‖2?w=?(wTXTXw?yTXw?wTXTy+yTy)?w=XTXw+(XTX)Tw?(yTX)T?XTy=2XTXw?2XTy
其中XTXw=(XTX)TwXTXw=(XTX)Tw.令求導(dǎo)結(jié)果等于0，即可以求導(dǎo)問題的最小值。
2XTXw?2XTy=0w=(XTX)?1XTy
2XTXw?2XTy=0w=(XTX)?1XTy
??再利用最小二乘法的矩陣形式對(duì)前面的例子進(jìn)行求解，用二階多項(xiàng)式擬合即兩個(gè)點(diǎn)(2,3),(5,8)(2,3),(5,8)。
表示輸入矩陣 XX和標(biāo)簽向量yy
X=[1125425]y=[38]T
X=[1241525]y=[38]T
計(jì)算XTXXTX
XTX=???272972913329133641???
XTX=[272972913329133641]
矩陣求逆，再做矩陣乘法運(yùn)算
但 XTXXTX不可逆，故無(wú)唯一解。
??關(guān)于矩陣的逆是否存在，可以通過(guò)判斷矩陣的行列式是否為0（det(A)=?0det(A)=?0 來(lái)判斷，也可以通過(guò)初等行變換，觀察矩陣的行向量是否線性相關(guān)，在這個(gè)例子下，矩陣不可逆，故有無(wú)窮多解。但如果新增一個(gè)點(diǎn)(4,7)(4,7)，則就可以解了。
??其實(shí)這和數(shù)據(jù)集的點(diǎn)數(shù)和選擇的階數(shù)有關(guān)，如果點(diǎn)數(shù)小于階數(shù)則會(huì)出現(xiàn)無(wú)窮解的情況，如果點(diǎn)數(shù)等于階數(shù)，那么剛好有解可以完全擬合所有數(shù)據(jù)點(diǎn)，如果點(diǎn)數(shù)大于階數(shù)，則會(huì)求的近似解。
??那么對(duì)于點(diǎn)數(shù)小于階數(shù)的情況，如何求解？在python的多項(xiàng)式擬合函數(shù)中是可以擬合的，而且效果不錯(cuò)，具體算法不是很了解，可以想辦法參考python的ployfit()函數(shù)的實(shí)現(xiàn)。
4. 擬合階數(shù)的選擇
?? 在前面的推導(dǎo)中，多項(xiàng)式的階數(shù)被固定了，那么實(shí)際場(chǎng)景下應(yīng)該如何選擇合適的階數(shù)MM呢？
一般會(huì)選擇階數(shù)MM小于點(diǎn)數(shù)NN
把訓(xùn)練數(shù)據(jù)分為訓(xùn)練集合驗(yàn)證集，在訓(xùn)練集上，同時(shí)用不同的MM值訓(xùn)練多個(gè)模型，然后選擇在驗(yàn)證集誤差最小的階數(shù)script type="math/tex" id="MathJax-Element-5573"M/script
numpy的allclose方法，比較兩個(gè)array是不是每一元素都相等，默認(rèn)在1e-05的誤差范圍內(nèi)。
使用如圖：
源碼如下：
python做科學(xué)計(jì)算的特點(diǎn)：1. 科學(xué)庫(kù)很全。（推薦學(xué)習(xí)：Python視頻教程）
科學(xué)庫(kù)：numpy，scipy。作圖：matplotpb。并行：mpi4py。調(diào)試：pdb。
2. 效率高。
如果你能學(xué)好numpy（array特性，f2py），那么你代碼執(zhí)行效率不會(huì)比f(wàn)ortran，C差太多。但如果你用不好array，那樣寫出來(lái)的程序效率就只能呵呵了。所以入門后，請(qǐng)一定花足夠多的時(shí)間去了解numpy的array類。
3. 易于調(diào)試。
pdb是我見過(guò)最好的調(diào)試工具，沒有之一。直接在程序斷點(diǎn)處給你一個(gè)截面，這只有文本解釋語(yǔ)言才能辦到。毫不夸張的說(shuō)，你用python開發(fā)程序只要fortran的1/10時(shí)間。
4. 其他。
它豐富而且統(tǒng)一，不像C++的庫(kù)那么雜（好比pnux的各種發(fā)行版），python學(xué)好numpy就可以做科學(xué)計(jì)算了。python的第三方庫(kù)很全，但是不雜。python基于類的語(yǔ)言特性讓它比起fortran等更加容易規(guī)?；_發(fā)。
數(shù)值分析中，龍格－庫(kù)塔法（Runge-Kutta methods）是用于非線性常微分方程的解的重要的一類隱式或顯式迭代法。這些技術(shù)由數(shù)學(xué)家卡爾·龍格和馬丁·威爾海姆·庫(kù)塔于1900年左右發(fā)明。
龍格-庫(kù)塔(Runge-Kutta)方法是一種在工程上應(yīng)用廣泛的高精度單步算法，其中包括著名的歐拉法，用于數(shù)值求解微分方程。由于此算法精度高，采取措施對(duì)誤差進(jìn)行抑制，所以其實(shí)現(xiàn)原理也較復(fù)雜。
高斯積分是在概率論和連續(xù)傅里葉變換等的統(tǒng)一化等計(jì)算中有廣泛的應(yīng)用。在誤差函數(shù)的定義中它也出現(xiàn)。雖然誤差函數(shù)沒有初等函數(shù)，但是高斯積分可以通過(guò)微積分學(xué)的手段解析求解。高斯積分（Gaussian integral），有時(shí)也被稱為概率積分，是高斯函數(shù)的積分。它是依德國(guó)數(shù)學(xué)家兼物理學(xué)家卡爾·弗里德里?！じ咚怪帐纤?。
洛倫茨吸引子及其導(dǎo)出的方程組是由愛德華·諾頓·洛倫茨于1963年發(fā)表，最初是發(fā)表在《大氣科學(xué)雜志》（Journal of the Atmospheric Sciences）雜志的論文《Deterministic Nonperiodic Flow》中提出的，是由大氣方程中出現(xiàn)的對(duì)流卷方程簡(jiǎn)化得到的。
這一洛倫茨模型不只對(duì)非線性數(shù)學(xué)有重要性，對(duì)于氣候和天氣預(yù)報(bào)來(lái)說(shuō)也有著重要的含義。行星和恒星大氣可能會(huì)表現(xiàn)出多種不同的準(zhǔn)周期狀態(tài)，這些準(zhǔn)周期狀態(tài)雖然是完全確定的，但卻容易發(fā)生突變，看起來(lái)似乎是隨機(jī)變化的，而模型對(duì)此現(xiàn)象有明確的表述。
更多Python相關(guān)技術(shù)文章，請(qǐng)?jiān)L問Python教程欄目進(jìn)行學(xué)習(xí)！以上就是小編分享的關(guān)于python能做什么科學(xué)計(jì)算的詳細(xì)內(nèi)容希望對(duì)大家有所幫助，更多有關(guān)python教程請(qǐng)關(guān)注環(huán)球青藤其它相關(guān)文章！
Python math 庫(kù)提供許多對(duì)浮點(diǎn)數(shù)的數(shù)學(xué)運(yùn)算函數(shù)，math模塊不支持復(fù)數(shù)運(yùn)算，若需計(jì)算復(fù)數(shù)，可使用cmath模塊（本文不贅述）。
使用dir函數(shù)，查看math庫(kù)中包含的所有內(nèi)容：
1） math.pi????# 圓周率π
2） math.e????#自然對(duì)數(shù)底數(shù)
3） math.inf? ? #正無(wú)窮大∞，-math.inf? ? #負(fù)無(wú)窮大-∞
4） math.nan? ? #非浮點(diǎn)數(shù)標(biāo)記，NaN（not a number）
1） math.fabs(x)? ? #表示X值的絕對(duì)值
2） math.fmod(x,y)? ? #表示x/y的余數(shù)，結(jié)果為浮點(diǎn)數(shù)
3） math.fsum([x,y,z])? ? #對(duì)括號(hào)內(nèi)每個(gè)元素求和，其值為浮點(diǎn)數(shù)
4） math.ceil(x)? ? #向上取整，返回不小于x的最小整數(shù)
5）math.floor(x)? ? #向下取整，返回不大于x的最大整數(shù)
6） math.factorial(x)? ? #表示X的階乘，其中X值必須為整型，否則報(bào)錯(cuò)
7） math.gcd(a,b)? ? #表示a,b的最大公約數(shù)
8）? math.frexp(x)? ? ? #x = i *2^j，返回（i，j）
9） math.ldexp(x,i)? ? #返回x*2^i的運(yùn)算值，為math.frexp(x)函數(shù)的反運(yùn)算
10） math.modf(x)? ? #表示x的小數(shù)和整數(shù)部分
11） math.trunc(x)? ? #表示x值的整數(shù)部分
12） math.copysign(x,y)? ? #表示用數(shù)值y的正負(fù)號(hào)，替換x值的正負(fù)號(hào)
13） math.isclose(a,b,rel_tol =x,abs_tol = y)? ? #表示a，b的相似性，真值返回True，否則False；rel_tol是相對(duì)公差：表示a，b之間允許的最大差值，abs_tol是最小絕對(duì)公差，對(duì)比較接近于0有用，abs_tol必須至少為0。
14） math.isfinite(x)? ? #表示當(dāng)x不為無(wú)窮大時(shí)，返回True，否則返回False
15） math.isinf(x)? ? #當(dāng)x為±∞時(shí)，返回True，否則返回False
16） math.isnan(x)? ? #當(dāng)x是NaN，返回True，否則返回False
1） math.pow(x,y)? ? #表示x的y次冪
2） math.exp(x)? ? #表示e的x次冪
3） math.expm1(x)? ? #表示e的x次冪減1
4） math.sqrt(x)? ? #表示x的平方根
5） math.log（x,base）? ? #表示x的對(duì)數(shù)值，僅輸入x值時(shí)，表示ln(x)函數(shù)
6） math.log1p(x)? ? #表示1+x的自然對(duì)數(shù)值
7） math.log2(x)? ? #表示以2為底的x對(duì)數(shù)值
8） math.log10(x)? ? #表示以10為底的x的對(duì)數(shù)值
1） math.degrees(x)? ? #表示弧度值轉(zhuǎn)角度值
2） math.radians(x)? ? #表示角度值轉(zhuǎn)弧度值
3） math.hypot(x,y)? ? #表示(x,y)坐標(biāo)到原點(diǎn)(0,0)的距離
4） math.sin(x)? ? #表示x的正弦函數(shù)值
5） math.cos(x)? ? #表示x的余弦函數(shù)值
6） math.tan(x)? ? #表示x的正切函數(shù)值
7）math.asin(x)? ? #表示x的反正弦函數(shù)值
8）?math.acos(x)? ? #表示x的反余弦函數(shù)值
9）?math.atan(x)? ? #表示x的反正切函數(shù)值
10） math.atan2(y,x)? ? #表示y/x的反正切函數(shù)值
11） math.sinh(x)? ? #表示x的雙曲正弦函數(shù)值
12） math.cosh(x)? ? #表示x的雙曲余弦函數(shù)值
13） math.tanh(x)? ? #表示x的雙曲正切函數(shù)值
14） math.asinh(x)? ? #表示x的反雙曲正弦函數(shù)值
15） math.acosh(x)? ? #表示x的反雙曲余弦函數(shù)值
16） math.atanh(x)? ? #表示x的反雙曲正切函數(shù)值
1）math.erf(x)? ? #高斯誤差函數(shù)
2） math.erfc(x)? ? #余補(bǔ)高斯誤差函數(shù)
3） math.gamma(x)? ? #伽馬函數(shù)（歐拉第二積分函數(shù)）
4） math.lgamma(x)? ? #伽馬函數(shù)的自然對(duì)數(shù)
前文提到了神經(jīng)網(wǎng)絡(luò)中的Sigmoid函數(shù)，實(shí)際上在反向傳播中還會(huì)用到Sigmoid的導(dǎo)數(shù)，形式很簡(jiǎn)單: s(x)*(1-s(x))，但是我想把這個(gè)過(guò)程自己推導(dǎo)一次，順便復(fù)習(xí)一下導(dǎo)數(shù)和微分。
Derivative(導(dǎo)數(shù))和Differential(微分)
首先我畫了一張圖來(lái)說(shuō)明什么是導(dǎo)數(shù)和微分，本質(zhì)上就是在極限中以線性函數(shù)(直線)來(lái)表示非線性函數(shù)(曲線)。
紅色的線是第一條割線(從[x,f(x)]到[x+h, f(x+h)])，(f(x+h) - f(x))/h 就是割線的斜率，物理學(xué)上是一段時(shí)間內(nèi)的平均速度。
灰色的線是第二條割線，當(dāng)割線圍繞著[x, f(x)]為原點(diǎn)繼續(xù)順時(shí)針轉(zhuǎn)動(dòng)時(shí)，h會(huì)不斷變小，小到極限就變成了[x, f(x)]的切線。
藍(lán)色的線即這條切線，其斜率就是[x,f(x)]的 導(dǎo)數(shù) ，物理意義是當(dāng)前這一個(gè)點(diǎn)的瞬間速度。
當(dāng)h小到極限的時(shí)候dy(導(dǎo)數(shù)除以h)就是[x,f(x)]的微分。
割線斜率減去切線斜率即為誤差函數(shù)E(h)
Reciprocal Rule(倒數(shù)法則)
根據(jù)微積分中的倒數(shù)法則，如果g(x) = 1/f(x), 則有
這個(gè)簡(jiǎn)單公式也非常容易證明
再將極限表達(dá)式分拆一下
因?yàn)閒在x點(diǎn)的連續(xù)性第二個(gè)極限表達(dá)式的分母等于f(x)的平方
現(xiàn)在利用倒數(shù)法則把Sigmoid函數(shù)的導(dǎo)數(shù)推導(dǎo)一下，這次我們記Sigmoid函數(shù)為s(x)，它的倒置函數(shù)為f(x)
根據(jù)倒數(shù)法則從f(x)開始推導(dǎo)得出公式S1
Chain Rule(鏈?zhǔn)椒▌t)
根據(jù)鏈?zhǔn)椒▌t我們可以有關(guān)于冪指求導(dǎo)的推廣
即
于是可以得出f(x)導(dǎo)數(shù)的另一種表達(dá)式S2
最后我們把S2和S1放到一起來(lái)消元就可以得到Sigmoid的導(dǎo)數(shù)公式了
用Python來(lái)實(shí)現(xiàn)如下邏輯:
References:
1. Differential on wiki
2. Chain rule on wiki
3.? Derivatives of logarithmic and exponential functions
4. MIT open course -?Multivariable Calculus
5.? Mathematics Stack Exchange
重慶分公司
重慶分公司
