Python 數(shù)據(jù)分析與科學(xué)計算 課件 第10章 Scipy科學(xué)計算

Python數(shù)據(jù)分析與科學(xué)計算第10章Scipy科學(xué)計算10.1Scipy中的科學(xué)計算工具10.2Scipy中的優(yōu)化10.3Scipy中的稀疏矩陣處理10.4典型應(yīng)用

常用SciPy模塊模塊說明模塊說明cluster向量量化optimize優(yōu)化constants數(shù)學(xué)常量linalg線性代數(shù)integrate積分stats統(tǒng)計interpolate插值io數(shù)據(jù)輸入輸出SciPy科學(xué)計算SciPy是Python中科學(xué)計算的核心包，SciPy以Numpy模塊為基礎(chǔ)添加了科學(xué)計算的模塊，包括最優(yōu)化、線性代數(shù)、積分、擬合、特殊函數(shù)、快速傅里葉變換、信號處理、圖像處理、常微分方程求解器等。

10.1 Scipy中的科學(xué)計算工具1.線性方程組求解SciPy提供了線性代數(shù)計算函數(shù)，這些函數(shù)大多放在scipy.linalg模塊中，可分為基本求解方法、特征值問題、矩陣分解、矩陣函數(shù)、矩陣方程求解、特殊矩陣構(gòu)造等。例如：

求解線性方程使用scipy.linalg模塊中的solve()函數(shù)，該函數(shù)接收兩個參數(shù)：數(shù)組a和數(shù)組b。其中數(shù)組a表示系數(shù)，數(shù)組b表示等號右側(cè)值?！纠?0.1】求解下列線性方程組。importnumpyasnpfromscipyimportlinalga=np.array([[3,5,0],[1,-1,0],[0,3,1]])#系數(shù)矩陣b=np.array([5,3,-1])x=linalg.solve(a,b)#求解線性方程組print(x)10.1 Scipy中的科學(xué)計算工具該線性方程組用矩陣可表示為：

程序運行結(jié)果：[2.5-0.50.5]1.線性方程組求解

2.范數(shù)范數(shù)(norm)是數(shù)學(xué)中的一種基本概念。在泛函分析中，它定義在賦范線性空間中，并滿足一定的條件：①非負性；②齊次性；③三角不等式。常常被用來度量某個向量空間（或矩陣）中的每個向量的長度或大小。10.1 Scipy中的科學(xué)計算工具

Scipy提供了函數(shù)linalg.norm來計算向量或矩陣的范數(shù)。常用范數(shù)10.1 Scipy中的科學(xué)計算工具2.范數(shù)向量范數(shù)含義矩陣范數(shù)含義0范數(shù)向量中非零元素的個數(shù)1范數(shù)矩陣中每列元素和的最大值1范數(shù)向量元素絕對值之和-1范數(shù)矩陣中每列元素和的最小值2范數(shù)歐幾里得范數(shù)‘fro’范數(shù)矩陣中所有元素平方和的平方根∞范數(shù)所有向量元素絕對值中的最大值∞范數(shù)矩陣中每行元素和的最大值-∞范數(shù)所有向量元素絕對值中的最小值-∞范數(shù)矩陣中每行元素和的最小值p范數(shù)向量元素絕對值的p次方和的1/p次冪2范數(shù)矩陣的最大奇異值【例10.2】輸出矩陣中的常用范數(shù)。fromscipyimportlinalgimportnumpyasnpM=np.array([[5，6],[3,4]])#創(chuàng)建矩陣print('M矩陣的1范數(shù)為：\n',linalg.norm(M,1))print('M矩陣的2范數(shù)為：\n',linalg.norm(M,2))print('M矩陣的正無窮范數(shù)為：\n',linalg.norm(M,np.inf))print('M矩陣的fro范數(shù)為：\n',linalg.norm(M,'fro'))10.1 Scipy中的科學(xué)計算工具2.范數(shù)程序運行結(jié)果：M矩陣的1范數(shù)為：10.0M矩陣的2范數(shù)為：9.271109059321276M矩陣的正無窮范數(shù)為：11.0M矩陣的fro范數(shù)為：9.273618495495704

3.統(tǒng)計分布Scipy中的stats模塊包含了各種連續(xù)分布和離散分布模型。連續(xù)變量分布表示：概率密度函數(shù)PDF累積分布函數(shù)CDF離散變量分布表示：概率質(zhì)量函數(shù)PMF累積分布函數(shù)CDF10.1 Scipy中的科學(xué)計算工具

（1）伯努利分布伯努利分布是一種離散分布，又稱為兩點分布或0-1分布，它的隨機變量只取0或者1。進行一次事件試驗，該事件發(fā)生的概率為p，不發(fā)生的概率為1-p。

任何一個只有兩種結(jié)果的隨機現(xiàn)象都服從0-1分布。

例如拋硬幣觀察正反面，打槍上靶統(tǒng)計等。3.統(tǒng)計分布10.1 Scipy中的科學(xué)計算工具【例10.3】伯努利分布應(yīng)用：打槍上靶實驗。importnumpyasnpfromscipyimportstatsimportmatplotlib.pyplotaspltimportmatplotlibasmplmpl.rcParams['font.sans-serif']=['SimHei']plt.grid()x=np.arange(0,2,1)p=0.75#打槍上靶率pList=stats.bernoulli.pmf(x,p)print(pList)plt.plot(x,pList,marker='o',linestyle='None')plt.vlines(x,(0,0),pList)plt.xlabel('隨機變量：打槍上靶1次')plt.ylabel('概率')plt.title('伯努利分布：p=%.2f'%p)plt.show()程序運行結(jié)果：[0.250.75]10.1 Scipy中的科學(xué)計算工具

（2）二項分布又稱為n重伯努利分布。如果隨機變量序列Xn(n=1,2,…)中的隨機變量均服從參數(shù)為p的伯努利分布，那么隨機變量序列Xn就形成了參數(shù)為p的n重伯努利分布。

例如，假定重復(fù)拋擲一枚均勻硬幣n次，如果在第i次拋擲中出現(xiàn)正面，令Xi=1；如果出現(xiàn)反面，則令Xi=0。那么，隨機變量Xn(n=1,2,…)就形成了參數(shù)為1/2的n重伯努利試驗。10.1 Scipy中的科學(xué)計算工具3.統(tǒng)計分布【例10.4】二項分布應(yīng)用：打槍上靶實驗，打靶10次，命中8環(huán)的概率分布。importnumpyasnpfromscipyimportstatsimportmatplotlib.pyplotaspltimportmatplotlibasmplmpl.rcParams['font.sans-serif']=['SimHei']plt.grid()n=10#打搶次數(shù)p=0.75#8環(huán)命中率10.1 Scipy中的科學(xué)計算工具x=np.arange(0,n+1,1)pList=stats.binom.pmf(x,n,p)print(pList)plt.plot(x,pList,marker='o’,linestyle='None')plt.vlines(x,0,pList)plt.xlabel('隨機變量：打中八環(huán)次數(shù)')plt.ylabel('二項分布概率')plt.title('二項分布：n=%i,p=%.3f'%(n,p))plt.show()程序運行結(jié)果：[9.53674316e-072.86102295e-053.86238098e-043.08990479e-031.62220001e-025.83992004e-021.45998001e-012.50282288e-012.81567574e-011.87711716e-015.63135147e-02]10.1 Scipy中的科學(xué)計算工具

(3)正態(tài)分布正態(tài)分布也稱為高斯分布，是統(tǒng)計學(xué)中最重要的連續(xù)概率分布。若隨機變量X服從一個數(shù)學(xué)期望為μ、方差為δ^2的正態(tài)分布，記為N(μ，δ^2)。其概率密度函數(shù)為正態(tài)分布的期望值μ決定了其位置，其標準差δ決定了分布的幅度。

當μ=0、δ=1時的正態(tài)分布是標準正態(tài)分布。10.1 Scipy中的科學(xué)計算工具3.統(tǒng)計分布【例10.5】正態(tài)分布應(yīng)用。importnumpyasnpfromscipyimportstatsimportmatplotlib.pyplotaspltimportmatplotlibasmplmpl.rcParams['font.sans-serif']=['SimHei']mu=0#平均值sigma=1#標準差x=np.arange(-10,10,0.1)y=stats.norm.pdf(x,mu,sigma)plt.plot(x,y)plt.xticks(np.arange(-10,10,1))plt.xlabel('隨機變量：x')plt.ylabel('概率')plt.title('正態(tài)分布：$\mu$=%.1f,$\sigma^2$=%.1f'%(mu,sigma))plt.grid()plt.show()10.1 Scipy中的科學(xué)計算工具

4.積分egrate模塊是SciPy中負責(zé)處理積分操作的模塊。10.1 Scipy中的科學(xué)計算工具egrate模塊中的常用函數(shù)函數(shù)說明quad一重積分dblquad二重積分tplquad三重積分nquadn倍多重積分fixed_quad高斯積分【例10.6】使用Segrate實現(xiàn)[0,1]區(qū)間上高斯函數(shù)的積分。importegratefromnumpyimportexpf=lambdax:exp(-x**2)i=egrate.quad(f,0,1)print(i)10.1 Scipy中的科學(xué)計算工具程序運行結(jié)果：(0.7468241328124271,8.291413475940725e-15)【例10.7】使用Segrate計算二重積分。importegratefrommathimportsqrtf=lambdax,y:8*x*yg=lambdax:0h=lambday:sqrt(1-2*y**2)i=egrate.dblquad(f,0,0.5,g,h)print(i)10.1 Scipy中的科學(xué)計算工具程序運行結(jié)果：(0.375,1.2071790918144718e-14)5.插值10.1 Scipy中的科學(xué)計算工具插值是通過已知的離散點求未知數(shù)據(jù)的過程或方法。SciPy提供了erpolate模塊進行插值操作：fromscipyimportinterpolate（1）一維插值

在給定有限數(shù)據(jù)點集合的情況下，通過構(gòu)建一個函數(shù)來近似估計這些數(shù)據(jù)點之間的值。

interp1d()函數(shù)：一維數(shù)據(jù)的插值運算

該函數(shù)接收兩個參數(shù)x點和y點，返回值是可調(diào)用函數(shù)

該函數(shù)可以用新的x調(diào)用并返回相應(yīng)的y，y=f(x)?！纠?0.8】一維插值應(yīng)用。fromerpolateimportinterp1dimportnumpyasnpimportmatplotlib.pyplotaspltimportrandomnp.random.seed(200)x=np.random.randn(10)y=np.random.randn(10)f=interp1d(x,y,kind='cubicx_new=np.linspace(start=x.min(),stop=x.max(),num=40)y_new=f(x_new)plt.plot(x,y,'ro',x_new,y_new,'*b-')plt.show()10.1 Scipy中的科學(xué)計算工具5.插值10.1 Scipy中的科學(xué)計算工具（2）樣條插值在樣條插值中，將點與由多項式定義的分段函數(shù)(稱為樣條)擬合。

InterpolatedUnivariateSpline()：樣條插值函數(shù)

接受x和y作為參數(shù)并生成可調(diào)用的函數(shù)，該函數(shù)可通過新的x進行調(diào)用?！纠?0.9】樣條插值應(yīng)用。fromerpolateimportInterpolatedUnivariateSplineimportnumpyasnpimportmatplotlib.pyplotaspltimportrandomnp.random.seed(10000)x=np.linspace(-5,5,20)y=np.exp(-x**2)+0.01*np.random.randn(20)+1f=InterpolatedUnivariateSpline(x,y)x_new=np.linspace(-5,5,100)y_new=f(x_new)plt.plot(x,y,'r*',x_new,y_new,'b-',lw=2)plt.show()10.1 Scipy中的科學(xué)計算工具

1.方程求解及求極值optimize.root()函數(shù)：求解非線性方程的根函數(shù)參數(shù)：

fun：表示方程的函數(shù)。

x0：根的初始猜測。10.2 Scipy中的優(yōu)化10.2 Scipy中的優(yōu)化1.方程求解及求極值optimize.minimize()函數(shù)：求函數(shù)最小值函數(shù)參數(shù)：

fun：要優(yōu)化的函數(shù)

x0：初始猜測值

method：要使用的方法名稱，值可以是：'CG'，'BFGS'，'Newton-CG'，'L-BFGS-B'，'TNC'，'COBYLA'，'SLSQP'。【例10.10】使用root()函數(shù)求解方程

的根。fromscipy.optimizeimportrootfrommathimportcosdeffunc(x):returnx+cos(x)myroot=root(func,0)print(myroot.x)10.2 Scipy中的優(yōu)化程序運行結(jié)果：[-0.73908513]

importnumpyasnpfromscipyimportoptimizeimportmatplotlib.pyplotaspltdeff(x):returnx**2+3*x+4x=np.arange(-10,10,0.1)plt.plot(x,f(x))#繪制函數(shù)曲線xout=optimize.fmin_bfgs(f,0)xmin=xout[0]ymin=f(xmin)print('xmin:',xmin)print('ymin:',ymin)plt.scatter(xmin,ymin,s=15,c='r’)plt.show()10.2 Scipy中的優(yōu)化程序運行結(jié)果：Optimizationterminatedsuccessfully.Currentfunctionvalue:1.750000Iterations:2Functionevaluations:6Gradientevaluations:3xmin:-1.5000000141634802ymin:1.7510.2 Scipy中的優(yōu)化函數(shù)曲線及最小值2.數(shù)據(jù)擬合擬合是用一個連續(xù)函數(shù)（曲線）靠近給定的離散數(shù)據(jù)，使其與給定的數(shù)據(jù)相吻合。(1)多項式擬合用一個n階多項式描述點(x,y)的關(guān)系：10.2 Scipy中的優(yōu)化

多項式擬合的目的是找到一組系數(shù)a，使得擬合得到的曲線和真實數(shù)據(jù)點之間的距離最小。多項式擬合的系數(shù)可以使用Numpy模塊的np.ployfit()函數(shù)得到?！纠?0.12】多項式擬合。importnumpyasnpimportmatplotlib.pyplotaspltx=np.linspace(-5,5,100)y=4*x+15#待擬合函數(shù)y_noise=y+np.random.randn(100)*2coeff=np.polyfit(x,y_noise,1)plt.plot(x,y_noise,'*r',x,coeff[0]*x+coeff[1])10.2 Scipy中的優(yōu)化2.數(shù)據(jù)擬合擬合結(jié)果如圖10.7所示。10.2 Scipy中的優(yōu)化2.數(shù)據(jù)擬合(2)曲線擬合曲線擬合可分為沒有約束的曲線擬合和帶有約束條件的曲線擬合。Scipy的optimize模塊中提供了curve_fit()函數(shù)來擬合沒有約束條件的曲線，least_squares()函數(shù)用來擬合帶有約束條件的曲線?！纠?0.13】曲線擬合。importnumpyasnpfromscipy.optimizeimportcurve_fitimportmatplotlib.pyplotaspltnp.random.seed(10)deffun(x):measurement_error=np.random.randn(len(x))y=5*x**2+measurement_errorreturnydefobjective_fun(x,a):returna*x**210.2 Scipy中的優(yōu)化x=np.arange(0,10,1)*0.2y=fun(x)pout,m=curve_fit(objective_fun,x,y)a=poutplt.scatter(x,y,label='data')plt.plot(x,objective_fun(x,a),"--",color="r",label="objectivequadraticfunction")plt.legend()plt.show()【例10.13】曲線擬合。10.2 Scipy中的優(yōu)化曲線擬合結(jié)果

10.3 Scipy中的稀疏矩陣處理在矩陣中，若數(shù)據(jù)為零的元素數(shù)目遠遠多于非零元素的數(shù)目，并且非零元素分布沒有規(guī)律時，則稱該矩陣為稀疏矩陣。根據(jù)存儲方式不同，可以將稀疏矩陣分成以下幾類：bsr_matrix：分塊存儲，基于行coo_matrix：坐標形式存儲（COO）csr_matrix：基于行的壓縮存儲（CSR）csc_matrix：基于列的壓縮存儲（CSC）dia_matrix：對角線存儲dok_matrix：基于鍵值對的存儲lil_matrix：基于行的鏈表存儲【例10.14】使用coo_matrix創(chuàng)建稀疏矩陣。importscipy.sparseasspvalues=[3,2,1,0]rows=[0,1,2,3]cols=[1,3,2,0]A=sp.coo_matrix((values,(rows,cols)),shape=[4,4])A.toarrayprint(A)10.3 Scipy中的稀疏矩陣處理程序運行結(jié)果：

(0,1) 3(1,3) 2(2,2) 1(3,0) 0【例10.15】存儲稀疏矩陣。importnumpyasnpfromscipy.sparseimportcsr_matrixindptr=np.array([0,0,3,4,4,5])indices=np.array([0,2,3,2,3,])data=np.array([8,1,7,1,2])csr=csr_matrix((data,indices,indptr))csr.toarray()10.3 Scipy中的稀疏矩陣處理程序運行結(jié)果：array([[0,0,0,0],[8,0,1,7],[0,0,1,0],[0,0,0,0],[0,0,0,2]])【例10.16】使用data屬性查看存儲的數(shù)據(jù)（非零項）。importnumpyasnpfromscipy.sparseimportcsr_matrixarr=np.array([[0,0,0,0],[1,0,0,0],[2,0,2,0],[0,0,0,0]])print(csr_matrix(arr).data)10.3 Scipy中的稀疏矩陣處理程序程序運行結(jié)果：[122]【例10.17】使用count_nonzero()方法計數(shù)非零元素的總數(shù)。importnumpyasnpfromscipy.sparseimportcsr_matrixarr=np.array([[0,0,0,0],[1,0,0,0],[2,0,2,0],[0,0,0,0]])print(csr_matrix(arr).count_nonzero())10.3 Scipy中的稀疏矩陣處理程序運行結(jié)果：3【例10.18】使用eliminate_zeros()方法從矩陣中刪除零項。importnumpyasnpfromscipy.sparseimportcsr_matrixarr=np.array([[0,0,0,0],[1,0,0,0],[2,0,2,0],[0,0,0,0]])m=csr_matrix(arr)m.eliminate_zeros()print(m)10.3 Scipy中的稀疏矩陣處理程序運行結(jié)果：

(1,0) 1(2,0) 2(2,2) 2【例10.19】使用sum_duplicates()方法消除重復(fù)項。importnumpyasnpfromscipy.sparseimportcsr_matrixarr=np.array([[0,0,0,0],[1,0,0,0],[2,0,2,0],[0,0,0,0]])mat=csr_matrix(arr)mat.sum_duplicates()print(mat)10.3 Scipy中的稀疏矩陣處理程序運行結(jié)果：

(1,0) 1(2,0) 2(2,2) 22.稀疏矩陣的運算Scipy.sparce提供了許多稀疏矩陣的運算，包括矩陣加法、減法、乘法以及矩陣求逆等。10.3 Scipy中的稀疏矩陣處理【例10.20】稀疏矩陣加法、減法和乘法運算。importnumpyasnpfromscipyimportsparseindptr=np.array([0,0,3,4,4,5])indices=np.array([0,2,3,2,3])data=np.array([8,1,7,1,2])csrA=sparse.csr_matrix((data,indices,indptr))print("矩陣A：\n",csrA.toarray())indptr1=np.array([0,1,3,3,4,5])indices1=np.array([0,2,3,2,3])10.3 Scipy中的稀疏矩陣處理data1=np.array([6,1,2,1,9])csrB=sparse.csr_matrix((data1,indices1,indptr1))print("矩陣B：\n",csrB.toarray())print("相加結(jié)果：\n",(csrB+csrA).toarray())print("相減結(jié)果：\n",(csrB-csrA).toarray())print("相乘結(jié)果：\n",sparse.csc_matrix.multiply(csrA,csrB).toarray())10.3 Scipy中的稀疏矩陣處理程序運行結(jié)果：矩陣A：

[[0000][8017][0010][0000][0002]]矩陣B：

[[6000][0012][0000][0010][0009]]相加結(jié)果：

[[6000][8029][0010][0010][00