python機器學習-乳腺癌細胞挖掘（博主親自錄制視頻）
機器學習，統計項目聯系：QQ：231469242
# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pylab as plt
from scipy import stats
x = [3.5, 2.5, 4.0, 3.8, 2.8, 1.9, 3.2, 3.7, 2.7, 3.3] #高中平均成績
y = [3.3, 2.2, 3.5, 2.7, 3.5, 2.0, 3.1, 3.4, 1.9, 3.7] #大學平均成績
#linregress(x,y)線性回歸函數
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
'''
Out[37]: LinregressResult(slope=0.70389344262295073, intercept=0.71977459016393475,
rvalue=0.68345387256609358, pvalue=0.029341978126562161, stderr=0.26581031503816904)
'''
slope = round(slope,3)
intercept = round(intercept,3)
print slope, intercept
#繪圖用
def f(x, a, b):
return a + b*x
xdata = np.linspace(1, 5, 20)
plt.grid(True)
plt.xlabel('x axis')
plt.ylabel('y axis')
plt.text(2.5, 4.0, r'$y = ' + str(intercept) + ' + ' + str(slope) +'*x$', fontsize=18)
plt.plot(xdata, f(xdata, intercept,slope), 'b', linewidth=1)
plt.plot(x,y,'ro')
plt.show()
作者Toby qq:231469242
http://www.cnblogs.com/NanShan2016/p/5493429.html
最小二乘法本身實現起來也是不難的，就如我們上面所說的不斷調整參數，然后令誤差函數Err不斷減小就行了。
With the advent of cheap computing power, statistical modeling has been a booming
field. This has also affected classical statistical analysis, as most problems can be
viewed from two perspectives: one can either make a statistical hypothesis, and
verify or falsify that hypothesis; or one can make a statistical model, and analyze
the significance of the model parameters.
Let me use a classical t-test as an example.
StatisticalModeling
Expressed as a statistical model, we assume that the difference between the first
and the second race is simply a constant value. (The null hypothesis would be that
this value is equal to zero.) This model has one parameter: the constant value. We
can find this parameter, as well as its confidence interval and a lot of additional
information, with the following Python code
# -*- coding: utf-8 -*-
'''
The command random.seed(123) initializes the random number generator with
the number 123, which ensures that two consecutive runs of this code produce the
same result, corresponding to the numbers given above.
'''
import numpy as np
from scipy import stats
# Generate the data
np.random.seed(123)
race_1 = np.round(np.random.randn(20)*10+90)
race_2 = np.round(np.random.randn(20)*10+85)
# t-test
(t, pVal) = stats.ttest_rel (race_1, race_2)
# Show the result
print('The probability that the two distributions '
'are equal is {0:5.3f} .'.format(pVal))
import pandas as pd
import statsmodels.formula.api as sm
np.random.seed(123)
df = pd.DataFrame({'Race1': race_1, 'Race2':race_2})
result = sm.ols(formula='I(Race2-Race1) ~ 1', data=df).fit()
print(result.summary())
statsmodels是一個包含統計模型、統計測試和統計數據挖掘python模塊。對每一個模型都會生成一個對應的統計結果。統計結果會和現有的統計包進行對比來保證其正確性。
1829年，高斯提供了最小二乘法的優化效果強於其他方法的證明，因此被稱為高斯-馬爾可夫定理。（來自於wikipedia）
最小二乘法(Least Squares Method,簡記為LSE)是一個比較古老的方法，源於天文學和測地學上的應用需要。在早期數理統計方法的發展中，這兩門科學起了很大的作用。丹麥統計學家霍爾把它們稱為“數理統計學的母親”。此后近三百年來，它廣泛應用於科學實驗與工程技術中。美國統計史學家斯蒂格勒( S. M. Stigler)指出, 最小二乘方法是19世紀數理統計學的壓倒一切的主題。1815年時,這方法已成為法國、意大利和普魯士在天文和測地學中的標准工具,到1825年時已在英國普遍使用。
追溯到1801年，意大利天文學家朱賽普·皮亞齊發現了第一顆小行星谷神星。經過40天的跟蹤觀測后，由於谷神星運行至太陽背后，使得皮亞齊失去了谷神星的位置。隨后全世界的科學家利用皮亞齊的觀測數據開始尋找谷神星，但是根據大多數人計算的結果來尋找谷神星都沒有結果。時年24歲的高斯也計算了谷神星的軌道。奧地利天文學家海因里希·奧爾伯斯根據高斯計算出來的軌道重新發現了谷神星。高斯於其1809年的著作《關於繞日行星運動的理論》中。在此書中聲稱他自1799年以來就使用最小二乘方法,由此爆發了一場與勒讓德的優先權之爭。
近代學者經過對原始文獻的研究,認為兩人可能是獨立發明了這個方法,但首先見於書面形式的，以勒讓德為早。然而,現今教科書和著作中,多把這個發明權歸功於高斯。其原因,除了高斯有更大的名氣外,主要可能是因為其正態誤差理論對這個方法的重要意義。勒讓德在其著作中,對最小二乘方法的優點有所闡述。然而,缺少誤差分析。我們不知道,使用這個方法引起的誤差如何,就需建立一種誤差分析理論。高斯於1823年在誤差e1 ,… , en獨立同分布的假定下,證明了最小二乘方法的一個最優性質: 在所有無偏的線性估計類中,最小二乘方法是其中方差最小的！在德國10馬克的鈔票上有高斯像，並配了一條正態曲線。在高斯眾多偉大的數學成就中挑選了這一條,亦可見這一成就對世界文明的影響。
現行的最小二乘法是勒讓德( A. M. Legendre)於1805年在其著作《計算慧星軌道的新方法》中提出的。它的主要思想就是選擇未知參數，使得理論值與觀測值之差的平方和達到最小：
我們現在看來會覺得這個方法似乎平淡無奇,甚至是理所當然的。這正說明了創造性思維之可貴和不易。從一些數學大家未能在這個問題上有所突破,可以看出當時這個問題之困難。歐拉、拉普拉斯在許多很困難的數學問題上有偉大的建樹,但在這個問題上未能成功。
在高斯發表其1809年著作之前，約在1780年左右,拉普拉斯已發現了概率論中的“中心極限定理”。根據這個定理，大量獨立的隨機變量之和,若每個變量在和中起的作用都比較小，則和的分布必接近於正態。測量誤差正具有這種性質。一般地說,隨機(而非系統)的測量誤差，是出自大量不顯著的來源的疊加。因此，中心極限定理給誤差的正態性提供了一種合理的理論解釋。這一點對高斯理論的圓滿化很有意義,因為高斯原來的假定(平均數天然合理)總難免給人一種不自然的感覺。
耐人尋味的是,無論是中心極限定理的發明者拉普拉斯,還是早就了解這一結果的高斯,都沒有從這個結果的啟示中去考察誤差分布問題。對前者而言,可能是出於思維定勢的束縛,這對拉普拉斯來說可算不幸,他因此失掉了把這個重要分布冠以自己名字的機會(正態分布這個形式最早是狄莫弗( De Moiv re) 1730年在研究二項概率的近似計算時得出的。以后也有其他學者使用過,但都沒有被冠以他們的名字。高斯之所以獲得這一殊榮,無疑是因為他把正態分布與誤差理論聯系了起來) 。
可以說,沒有高斯的正態誤差理論配合, 最小二乘方法的意義和重要性可能還不到其現今所具有的十分之一。最小二乘方法方法與高斯誤差理論的結合,是數理統計史上最重大的成就之一,其影響直到今日也尚未過時！由於本文是主要介紹最小二乘法與矩陣投影之間的關系，對於最小二乘和概率之間的關系，請參看靳志輝的《正態分布的前世今生》。
那么，投影矩陣與最小二乘二者有什么必然的聯系么，當我開始寫這篇文章的時候我也這樣問自己。先說說投影吧，這個想必大家都知道，高中的知識了。一個向量在另一個向量上的投影，實際上就是尋找在上離最近的點。
現在我們假設投影點是向量上的一點p，可以規定p=xa（x是某個數）。定義e=b-p，稱e 為誤差。因為e 與p 也就是a 垂直，所以有
，展開化簡得到：
，
我們發現：如果改變b，那么p相對應改變，然而改變a，p無變化。接下來，我們可以考慮更高維度的投影，三維空間的投影是怎么樣的呢，我們可以想象一個三維空間內的向量在該空間內的一個平面上的投影：
我們假設這個平面的基(basis)是a1, a2。那么矩陣A 的列空間就是該平面。假設一個不在該平面上的向量b 在該平面上的投影是p 。我們的任務就是找到合適的x，使得p=Ax 。這里有一個關鍵的地方：e 與該平面垂直，所以
。我們把上邊式子展開，得到
,
有了上面的背景知識，我們可以正式進入主題了，投影矩陣(projection matrix)：
這里我們最需要關注的是投影矩陣的兩個性質：
1）
;
2）
;
對於第一個，很容易理解，因為P本身就是個對稱陣。第二個，直觀的理解就是投影到a上后再投影一次，顯然投影並沒有改變，也就是二次投影還是其本身。
這個投影到底有什么用呢？從上面的分析中我們可以看出：投影矩陣P可以吧向量b投影成向量p！從線性代數的角度來說，Ax=b並不一定總有解，這在實際情況中會經常遇到（m > n）。所以我們就把b投影到向量p上，因為p在a1，a2的平面內，所以Ax =p是可以求解的。
好了，在此我們先暫別“投影”。下面，開始說一下最小二乘的故事吧：在實際應用中，線性回歸是經常用到的，我們可以在一張散列點圖中作一條直線（暫時用直線）來近似表述這些散列點的關系。比如：
設變量y 與t 成線性關系,即.現在已知m 個實驗點ai和bi ,求兩個未知參數C,D 。將代入得矛盾方程組
令
,
則可寫成Ax=b形式。從線性代數的角度來看，就是A的列向量的線性組合無法充滿整個列空間，也就是說Ax=b這個方程根本沒有解。從圖形上也很好理解：根本沒有一條直線同時經過所有藍色的點！所以為了選取最合適的x，讓該等式"盡量成立"，引入殘差平方和函數H：
這也就是最小二乘法的思想。我們知道，當x取最優值的時候，Ax恰好對應圖中線上橙色的點，而b則對應圖中藍色的點，e的值則應紅色的線長。
看到這里你有沒有和之前投影的那部分知識聯系在一起呢？最小二乘的思想是想如何選取參數x使得H最小。而從向量投影的角度來看這個問題，H就是向量e長度的平方，如何才能使e的長度最小呢？b和a1，a2都是固定的，當然是e垂直a1，a2平面的時候長度最小！換句話說：最小二乘法的解與矩陣投影時對變量求解的目標是一致的！
為了定量地給出與實驗數據之間線性關系的符合程度,可以用相關系數來衡量.它定義為
r也就是我們之前介紹的向量夾角。r 值越接近1, y與t 的線性關系越好.為正時,直線斜率為正,稱為正相關；r 為負時,直線斜率為負,稱為負相關.接近於0時,測量數據點分散或之間為非線性.不論測量數據好壞都能求出和,所以我們必須有一種判斷測量數據好壞的方法,用來判斷什么樣的測量數據不宜擬合,判斷的方法是時,測量數據是非線性的. r0稱為相關系數的起碼值,與測量次數n 有關。
最小二乘講到這里似乎已經說完了，但是有一個問題，那就是我們所利用的投影矩陣P這里我們假定
是可逆的，這種假定合理嗎？Strang在最后給我們作了解答：
If A has independent columns, then A'A is invertible
寫到這里，我想有必要總結一下，為什么最小二乘和投影矩陣要扯到一起，它們有什么聯系：最小二乘是用於數據擬合的一個很霸氣的方法，這個擬合的過程我們稱之為線性回歸。如果數據點不存在離群點（outliers），那么該方法總是會顯示其簡單粗暴的一面。我們可以把最小二乘的過程用矩陣的形式描述出來，然而，精妙之處就在於，這與我們的投影矩陣的故事不謀而合，所以，我們又可以借助於投影矩陣的公式，也就是
來加以解決。
最小二乘法是從誤差擬合角度對回歸模型進行參數估計或系統辨識,並在參數估計、系統辨識以及預測、預報等眾多領域中得到極為廣泛的應用。在數據擬合領域，最小二乘法及其各種變形的擬合方法包括:一元線性最小二乘法擬合、多元線性擬合、多項式擬合、非線性擬合。最小二乘法能將從實驗中得出的一大堆看上去雜亂無章的數據中找出一定規律,擬合成一條曲線來反映所給數據點總趨勢,以消除其局部波動。它為科研工作者提供了一種非常方便實效的數據處理方法。隨着現代電子計算機的普及與發展，這個占老的方法更加顯示出其強大的生命力。
想了解更多有關矩陣的內容，可以搜索《神奇的矩陣》。
http://www.cnblogs.com/NanShan2016/p/5493429.html參考
自從開始做畢設以來，發現自己無時無刻不在接觸最小二乘法。從求解線性透視圖中的消失點，m元n次函數的擬合，包括后來學到的神經網絡，其思想歸根結底全都是最小二乘法。
1-1 “多線→一點”視角與“多點→一線”視角
最小二乘法非常簡單，我把它分成兩種視角描述：
（1）已知多條近似交匯於同一個點的直線，想求解出一個近似交點：尋找到一個距離所有直線距離平方和最小的點，該點即最小二乘解；
（2）已知多個近似分布於同一直線上的點，想擬合出一個直線方程：設該直線方程為y=kx+b，調整參數k和b，使得所有點到該直線的距離平方之和最小，設此時滿足要求的k=k0，b=b0，則直線方程為y=k0x+b0。
1-2 思維拓展
這只是舉了兩個簡單的例子，其實在現實生活中我們可以利用最小二乘法解決更為復雜 的問題。比方說有一個未知系數的二元二次函數f(x,y)=w0x^2+w1y^2+w2xy+w3x+w4y+w5，這里w0~w5為未知的參數，為了 確定下來這些參數，將會給定一些樣本點(xi,yi,f(xi,yi))，然后通過調整這些參數，找到這樣一組w0~w5，使得這些所有的樣本點距離函數 f(x,y)的距離平方之和最小。至於具體用何種方法來調整這些參數呢？有一種非常普遍的方法叫“梯度下降法”，它可以保證每一步調整參數，都使得 f(x,y)朝比當前值更小的方向走，只要步長α選取合適，我們就可以達成這種目的。
而這里不得不提的就是神經網絡了。神經網絡其實就是不斷調整權值w和偏置b，來使 得cost函數最小，從這個意義上來講它還是屬於最小二乘法。更為可愛的一點是，神經網絡的調參用到的仍是梯度下降法，其中最常用的當屬隨機梯度下降法。 而后面偉大的bp算法，其實就是為了給梯度下降法做個鋪墊而已，bp算法的結果是cost函數對全部權值和全部偏置的偏導，而得知了這些偏導，對於各個權 值w和偏置b該走向何方就指明了方向。
因此，最小二乘法在某種程度上無異於機器學習中基礎中的基礎，且具有相當重要的地位。至於上面所說的“梯度下降法”以及“利用最小二乘法求解二元二次函數的w0~w5”，我將會在后面的博客中進行更加詳細的探討。
2 scipy庫中的leastsq函數
當然，最小二乘法本身實現起來也是不難的，就如我們上面所說的不斷調整參數，然后令誤差函數Err不斷減小就行了。我們將在下一次博客中詳細說明如何利用梯度下降法來完成這個目標。
而在本篇博客中，我們介紹一個scipy 庫中的函數，叫leastsq，它可以省去中間那些具體的求解步驟，只需要輸入一系列樣本點，給出待求函數的基本形狀（如我剛才所說，二元二次函數就是一 種形狀——f(x,y)=w0x^2+w1y^2+w2xy+w3x+w4y+w5，在形狀給定后，我們只需要求解相應的系數w0~w6），即可得到相應 的參數。至於中間到底是怎么求的，這一部分內容就像一個黑箱一樣。
2-1 函數形為y=kx+b
這一次我們給出函數形y=kx+b。這種情況下，待確定的參數只有兩個：k和b。
此時給出7個樣本點如下：
1 Xi=np.array([8.19,2.72,6.39,8.71,4.7,2.66,3.78]) 2 Yi=np.array([7.01,2.78,6.47,6.71,4.1,4.23,4.05])
則使用leastsq函數求解其擬合直線的代碼如下：
1 ###最小二乘法試驗### 2 import numpy as np 3 from scipy.optimize import leastsq 4 5 ###采樣點(Xi,Yi)### 6 Xi=np.array([8.19,2.72,6.39,8.71,4.7,2.66,3.78]) 7 Yi=np.array([7.01,2.78,6.47,6.71,4.1,4.23,4.05]) 8 9 ###需要擬合的函數func及誤差error### 10 def func(p,x): 11 k,b=p 12 return k*x+b 13 14 def error(p,x,y,s): 15 print s 16 return func(p,x)-y #x、y都是列表，故返回值也是個列表 17 18 #TEST 19 p0=[100,2] 20 #print( error(p0,Xi,Yi) ) 21 22 ###主函數從此開始### 23 s="Test the number of iteration" #試驗最小二乘法函數leastsq得調用幾次error函數才能找到使得均方誤差之和最小的k、b 24 Para=leastsq(error,p0,args=(Xi,Yi,s)) #把error函數中除了p以外的參數打包到args中 25 k,b=Para[0] 26 print"k=",k,'\n',"b=",b 27 28 ###繪圖，看擬合效果### 29 import matplotlib.pyplot as plt 30 31 plt.figure(figsize=(8,6)) 32 plt.scatter(Xi,Yi,color="red",label="Sample Point",linewidth=3) #畫樣本點 33 x=np.linspace(0,10,1000) 34 y=k*x+b 35 plt.plot(x,y,color="orange",label="Fitting Line",linewidth=2) #畫擬合直線 36 plt.legend() 37 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的求解結果，第二個元素我暫時也不知道是什么意思，先留下來。
其擬合效果圖如下：
2-2 函數形為y=ax^2+bx+c
這一次我們給出函數形y=ax^2+bx+c。這種情況下，待確定的參數有3個：a，b和c。
此時給出7個樣本點如下：
1 Xi=np.array([0,1,2,3,-1,-2,-3]) 2 Yi=np.array([-1.21,1.9,3.2,10.3,2.2,3.71,8.7])
這一次的代碼與2-1差不多，除了把待求參數再增加一個，換了一下訓練樣本，換了一下func中給出的函數形，幾乎沒有任何變化。
1 ###最小二乘法試驗### 2 import numpy as np 3 from scipy.optimize import leastsq 4 5 ###采樣點(Xi,Yi)### 6 Xi=np.array([0,1,2,3,-1,-2,-3]) 7 Yi=np.array([-1.21,1.9,3.2,10.3,2.2,3.71,8.7]) 8 9 ###需要擬合的函數func及誤差error### 10 def func(p,x): 11 a,b,c=p 12 return a*x**2+b*x+c 13 14 def error(p,x,y,s): 15 print s 16 return func(p,x)-y #x、y都是列表，故返回值也是個列表 17 18 #TEST 19 p0=[5,2,10] 20 #print( error(p0,Xi,Yi) ) 21 22 ###主函數從此開始### 23 s="Test the number of iteration" #試驗最小二乘法函數leastsq得調用幾次error函數才能找到使得均方誤差之和最小的a~c 24 Para=leastsq(error,p0,args=(Xi,Yi,s)) #把error函數中除了p以外的參數打包到args中 25 a,b,c=Para[0] 26 print"a=",a,'\n',"b=",b,"c=",c 27 28 ###繪圖，看擬合效果### 29 import matplotlib.pyplot as plt 30 31 plt.figure(figsize=(8,6)) 32 plt.scatter(Xi,Yi,color="red",label="Sample Point",linewidth=3) #畫樣本點 33 x=np.linspace(-5,5,1000) 34 y=a*x**2+b*x+c 35 plt.plot(x,y,color="orange",label="Fitting Curve",linewidth=2) #畫擬合曲線 36 plt.legend() 37 plt.show()
不過我們發現，它依舊能夠非常順利地解出待求的三個參數。其擬合情況如圖所示：
2-3 leastsq擬合y=kx+b可視化
本部分內容是建立在2-1代碼的基礎上，用Mayavi繪3D圖，以簡單地說明最小二乘法到底是怎么一回事。該部分知識用到了mgrid函數，具體是如何實施的請移步《Python閑談（一）mgrid慢放》。
step 1：創建一個k矩陣和b矩陣。在mgrid擴展后，有：
（1）k=[k1,k2,k3,...,kn]
mgrid(k)（朝右擴展）=
[k1,k1,k1,...,k1]
[k2,k2,k2,...,k2]
[k3,k3,k3,...,k3]
...
[kn,kn,kn,...,kn]
（2）b=[b1,b2,b3,...,bn]
mgrid(b)（朝下擴展）=
[b1,b2,b3,...,bn]
[b1,b2,b3,...,bn]
[b1,b2,b3,...,bn]
...
[b1,b2,b3,...,bn]
其中k矩陣和b矩陣等大（皆為n維向量,或者說1*n的矩陣），且這兩個矩陣里面的元素都非常密集。舉個例子以說明什么叫矩陣中的元素很密集：a是個矩陣，假設aij 為a矩陣中第i行第j列元素，則aij 和 a{i+1}j 的差值很小，aij 和 ai{j+1} 的差值也很小。也就是同一行或者同一列中相鄰的兩個元素的值非常接近。為什么要讓矩陣元素如此密集呢？因為我們的根本目的是用“密集的離散”來逼近“連續”，這里的思想就像微積分一樣。
而放在這里，就是ku和k{u+1}很接近，bv和b{v+1}也很接近。
step 2：令k矩陣和b矩陣中的元素按照其位置一一對應。對應后的結果為：
Combine_kb=
[(k1,b1),(k1,b2),(k1,b3)...,(k1,bn)]
[(k2,b1),(k2,b2),(k2,b3)...,(k2,bn)]
[(k3,b1),(k3,b2),(k3,b3)...,(k3,bn)]
...
[(kn,b1),(kn,b2),(kn,b3)...,(kn,bn)]
step 3：對矩陣中每一個(ku,bv)，我們分別求出該種情況下每一個訓練樣本點的誤差平方之和，即有：
Err{(ku,bv)}=∑{i=1~m}（(yi-（ku*xi+bv）)**2）
其中m為給定的訓練樣本點的個數。例如在這里：
1 Xi=np.array([8.19,2.72,6.39,8.71,4.7,2.66,3.78]) 2 Yi=np.array([7.01,2.78,6.47,6.71,4.1,4.23,4.05])
則有m=7。
什么意思呢？舉個例子，當i=1的時候，這個時候把(x1,y1)（= （8.19,7.01））代入（(yi-（ku*xi+bv）)**2）式子里面，由於此時已經鎖定(ku,bv)，因此式中所有的數都是常數，我們可以 解出一個常數（(y1-（ku*x1+bv）)**2）。然后依次令i=2,3,4,...,7,可以分別求解出一個（(yi- （ku*xi+bv）)**2）值來，這7個（(yi-（ku*xi+bv）)**2）值加起來即Err{(ku,bv)}。
注意了，最終我們算出的那個Err{(ku,bv)}將會存放到ku、bv對應的那個位置，比方說u=3,v=2：
mgrid(k)=
[k1,k1,k1,...,k1]
[k2,k2,k2,...,k2]
[k3,k3,k3,...,k3]
...
[kn,kn,kn,...,kn]
mgrid(b)=
[b1,b2,b3,...,bn]
[b1,b2,b3,...,bn]
[b1,b2,b3,...,bn]
...
[b1,b2,b3,...,bn]
則剛才算出來的Err{(k3,b2)}應該放在這個位置：
Err=
[Err11,Err12,Err13,...,Err1n]
[Err21,Err22,Err23,...,Err2n]
[Err31,Err32,Err33,...,Err3n]
...
[Errn1,Errn2,Errn3,...,Errnn]
如此這般對於每一對(ku,bv)都這樣算，則上方的Err矩陣中每一個元素的值都可以算出來；將計算出的結果正確地放在Err矩陣中對應位置，即得到Err矩陣。
step 4：繪制曲面。
截至目前我們已經得到了兩個重要矩陣Combine_kb和Err，其中Combine_kb提供點的x、y軸坐標，Err矩陣提供點的z軸坐標。
Combine_kb=
[(k1,b1),(k1,b2),(k1,b3)...,(k1,bn)]
[(k2,b1),(k2,b2),(k2,b3)...,(k2,bn)]
[(k3,b1),(k3,b2),(k3,b3)...,(k3,bn)]
...
[(kn,b1),(kn,b2),(kn,b3)...,(kn,bn)]
Err=
[Err11,Err12,Err13,...,Err1n]
[Err21,Err22,Err23,...,Err2n]
[Err31,Err32,Err33,...,Err3n]
...
[Errn1,Errn2,Errn3,...,Errnn]
我們再將這兩個矩陣合並一下得到Combine_kbErr矩陣：
Combine_kbErr=
[(k1,b1,Err11),(k1,b2,Err12),(k1,b3,Err13)...,(k1,bn,Err1n)]
[(k2,b1,Err21),(k2,b2,Err22),(k2,b3,Err23)...,(k2,bn,Err2n)]
[(k3,b1,Err31),(k3,b2,Err32),(k3,b3,Err33)...,(k3,bn,Err3n)]
...
[(kn,b1,Errn1),(kn,b2,Errn2),(kn,b3,Errn3)...,(kn,bn,Errnn)]
在三維空間直角坐標系下繪制出Combine_kbErr中的每一個點，然后將這些點與其各自相鄰的點連起來，則得到我們想要的Err(k,b)函數曲面。
step 5：本部分代碼如下：
1 """part 2""" 2 ###定義一個函數，用於計算在k、b已知時∑((yi-(k*xi+b))**2)### 3 def S(k,b): 4 ErrorArray=np.zeros(k.shape) #k的shape事實上同時也是b的shape 5 for x,y in zip(Xi,Yi): #zip(Xi,Yi)=[(8.19,7.01),(2.72,2.78),...,(3.78,4.05)] 6 ErrorArray+=(y-(k*x+b))**2 7 return ErrorArray 8 9 ###繪制ErrorArray+最低點### 10 from enthought.mayavi import mlab 11 12 #畫整個Error曲面 13 k,b=np.mgrid[k0-1:k0+1:10j,b0-1:b0+1:10j] 14 Err=S(k,b) 15 face=mlab.surf(k,b,Err/500.0,warp_scale=1) 16 mlab.axes(xlabel='k',ylabel='b',zlabel='Error') 17 mlab.outline(face) 18 19 #畫最低點（即k，b所在處） 20 MinErr=S(k0,b0) 21 mlab.points3d(k0,b0,MinErr/500.0,scale_factor=0.1,color=(0.5,0.5,0.5)) #scale_factor用來指定點的大小 22 mlab.show()
對要點說明如下：
1、為了讓最小二乘法求解的結果出現在繪制曲面的范圍內，我們以最終leastsq求得的k0、b0為中心創建k向量和b向量。
2、傳入S函數的是k向量和b向量mgrid后的結果。
3、S函數中的ErrorArray+=(y-(k*x+b))**2 操作里，k、b皆為矩陣（是k、b向量mgrid后的結果），而x、y皆為常數，故這里的操作實際上是對矩陣的操作。這個ErrorArray就是上面我說的Err矩陣。
4、在繪圖時之所以對Err除以500，是因為Err和k、b的差距不是一般的大，直接繪圖會導致什么都看不出來。舉一個最簡單的例子就是比如我們要畫個二維直角坐標系下的圖，x的取值范圍是0~1，y的取值范圍是0~1000，而兩個坐標軸卻都按一個單位△x=△y=0.1來畫，想想看結果會成什么樣子？
這里也是同樣的道理，於是得給Err除以一個大數才能讓圖像正常顯示。
其實matplotlib畫三維坐標系下的圖會幫你調整到合適，只有Mayavi才會出現這種情況，反正注意一下比例問題就好了。
5、該程序除過繪制Err曲面外，還把(k0,b0)也畫出來了，見灰色小球。
step 6：整個程序的全部代碼如下，其中part1與2-1的代碼是完全一樣的。
1 ###【最小二乘法試驗】### 2 import numpy as np 3 from scipy.optimize import leastsq 4 5 ###采樣點(Xi,Yi)### 6 Xi=np.array([8.19,2.72,6.39,8.71,4.7,2.66,3.78]) 7 Yi=np.array([7.01,2.78,6.47,6.71,4.1,4.23,4.05]) 8 9 """part 1""" 10 ###需要擬合的函數func及誤差error### 11 def func(p,x): 12 k,b=p 13 return k*x+b 14 15 def error(p,x,y): 16 return func(p,x)-y #x、y都是列表，故返回值也是個列表 17 18 p0=[1,2] 19 20 ###最小二乘法求k0、b0### 21 Para=leastsq(error,p0,args=(Xi,Yi)) #把error函數中除了p以外的參數打包到args中 22 k0,b0=Para[0] 23 print"k0=",k0,'\n',"b0=",b0 24 25 """part 2""" 26 ###定義一個函數，用於計算在k、b已知時，∑((yi-(k*xi+b))**2)### 27 def S(k,b): 28 ErrorArray=np.zeros(k.shape) #k的shape事實上同時也是b的shape 29 for x,y in zip(Xi,Yi): #zip(Xi,Yi)=[(8.19,7.01),(2.72,2.78),...,(3.78,4.05)] 30 ErrorArray+=(y-(k*x+b))**2 31 return ErrorArray 32 33 ###繪制ErrorArray+最低點### 34 from enthought.mayavi import mlab 35 36 #畫整個Error曲面 37 k,b=np.mgrid[k0-1:k0+1:10j,b0-1:b0+1:10j] 38 Err=S(k,b) 39 face=mlab.surf(k,b,Err/500.0,warp_scale=1) 40 mlab.axes(xlabel='k',ylabel='b',zlabel='Error') 41 mlab.outline(face) 42 43 #畫最低點（即k，b所在處） 44 MinErr=S(k0,b0) 45 mlab.points3d(k0,b0,MinErr/500.0,scale_factor=0.1,color=(0.5,0.5,0.5)) #scale_factor用來指定點的大小 46 mlab.show()
3 結語
本次博客給出了最小二乘法的Python實現方法，它用到了scipy庫中的leastsq函數。在上面我們給出了兩個實例，分別實現了對一元一次函數的擬合和一元二次函數的擬合，而事實上，對於函數並不一定得是一元函數，對於更多元的函數也同樣能夠利用最小二乘法完成擬合工作，不過隨着元和次的增加，待求參數也就越來越多了，比方說二元二次函數就有6個待求參數w0~w6。
然為了更好地理解神經網絡的訓練算法，並不建議直接使用leastsq函數完成對未知參數的求解，因此在以后的博客中我會詳細說明如何利用梯度下降法來求解誤差函數的最小值。
4 下面要寫的博客
1、梯度下降法（什么是梯度下降法，如何使用梯度下降法求一元二次函數最小值，如何使用梯度下降法求二元二次函數的最小值）
2、最小二乘法擬合二元二次函數（即求解w0~w5，相當於對梯度下降法的一個應用）
以上兩個內容我會放在同一個博客中。
https://study.163.com/provider/400000000398149/index.htm?share=2&shareId=400000000398149（ 歡迎關注博主主頁，學習python視頻資源，還有大量免費python經典文章）
