Python 数据分析与科学计算 课件 第10章 Scipy科学计算

Python数据分析与科学计算第10章Scipy科学计算10.1Scipy中的科学计算工具10.2Scipy中的优化10.3Scipy中的稀疏矩阵处理10.4典型应用

常用SciPy模块模块说明模块说明cluster向量量化optimize优化constants数学常量linalg线性代数integrate积分stats统计interpolate插值io数据输入输出SciPy科学计算SciPy是Python中科学计算的核心包，SciPy以Numpy模块为基础添加了科学计算的模块，包括最优化、线性代数、积分、拟合、特殊函数、快速傅里叶变换、信号处理、图像处理、常微分方程求解器等。

10.1 Scipy中的科学计算工具1.线性方程组求解SciPy提供了线性代数计算函数，这些函数大多放在scipy.linalg模块中，可分为基本求解方法、特征值问题、矩阵分解、矩阵函数、矩阵方程求解、特殊矩阵构造等。例如：

求解线性方程使用scipy.linalg模块中的solve()函数，该函数接收两个参数：数组a和数组b。其中数组a表示系数，数组b表示等号右侧值。【例10.1】求解下列线性方程组。importnumpyasnpfromscipyimportlinalga=np.array([[3,5,0],[1,-1,0],[0,3,1]])#系数矩阵b=np.array([5,3,-1])x=linalg.solve(a,b)#求解线性方程组print(x)10.1 Scipy中的科学计算工具该线性方程组用矩阵可表示为：

程序运行结果：[2.5-0.50.5]1.线性方程组求解

2.范数范数(norm)是数学中的一种基本概念。在泛函分析中，它定义在赋范线性空间中，并满足一定的条件：①非负性；②齐次性；③三角不等式。常常被用来度量某个向量空间（或矩阵）中的每个向量的长度或大小。10.1 Scipy中的科学计算工具

Scipy提供了函数linalg.norm来计算向量或矩阵的范数。常用范数10.1 Scipy中的科学计算工具2.范数向量范数含义矩阵范数含义0范数向量中非零元素的个数1范数矩阵中每列元素和的最大值1范数向量元素绝对值之和-1范数矩阵中每列元素和的最小值2范数欧几里得范数‘fro’范数矩阵中所有元素平方和的平方根∞范数所有向量元素绝对值中的最大值∞范数矩阵中每行元素和的最大值-∞范数所有向量元素绝对值中的最小值-∞范数矩阵中每行元素和的最小值p范数向量元素绝对值的p次方和的1/p次幂2范数矩阵的最大奇异值【例10.2】输出矩阵中的常用范数。fromscipyimportlinalgimportnumpyasnpM=np.array([[5，6],[3,4]])#创建矩阵print('M矩阵的1范数为：\n',linalg.norm(M,1))print('M矩阵的2范数为：\n',linalg.norm(M,2))print('M矩阵的正无穷范数为：\n',linalg.norm(M,np.inf))print('M矩阵的fro范数为：\n',linalg.norm(M,'fro'))10.1 Scipy中的科学计算工具2.范数程序运行结果：M矩阵的1范数为：10.0M矩阵的2范数为：9.271109059321276M矩阵的正无穷范数为：11.0M矩阵的fro范数为：9.273618495495704

3.统计分布Scipy中的stats模块包含了各种连续分布和离散分布模型。连续变量分布表示：概率密度函数PDF累积分布函数CDF离散变量分布表示：概率质量函数PMF累积分布函数CDF10.1 Scipy中的科学计算工具

（1）伯努利分布伯努利分布是一种离散分布，又称为两点分布或0-1分布，它的随机变量只取0或者1。进行一次事件试验，该事件发生的概率为p，不发生的概率为1-p。

任何一个只有两种结果的随机现象都服从0-1分布。

例如抛硬币观察正反面，打枪上靶统计等。3.统计分布10.1 Scipy中的科学计算工具【例10.3】伯努利分布应用：打枪上靶实验。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()程序运行结果：[0.250.75]10.1 Scipy中的科学计算工具

（2）二项分布又称为n重伯努利分布。如果随机变量序列Xn(n=1,2,…)中的随机变量均服从参数为p的伯努利分布，那么随机变量序列Xn就形成了参数为p的n重伯努利分布。

例如，假定重复抛掷一枚均匀硬币n次，如果在第i次抛掷中出现正面，令Xi=1；如果出现反面，则令Xi=0。那么，随机变量Xn(n=1,2,…)就形成了参数为1/2的n重伯努利试验。10.1 Scipy中的科学计算工具3.统计分布【例10.4】二项分布应用：打枪上靶实验，打靶10次，命中8环的概率分布。importnumpyasnpfromscipyimportstatsimportmatplotlib.pyplotaspltimportmatplotlibasmplmpl.rcParams['font.sans-serif']=['SimHei']plt.grid()n=10#打抢次数p=0.75#8环命中率10.1 Scipy中的科学计算工具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('随机变量：打中八环次数')plt.ylabel('二项分布概率')plt.title('二项分布：n=%i,p=%.3f'%(n,p))plt.show()程序运行结果：[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中的科学计算工具

(3)正态分布正态分布也称为高斯分布，是统计学中最重要的连续概率分布。若随机变量X服从一个数学期望为μ、方差为δ^2的正态分布，记为N(μ，δ^2)。其概率密度函数为正态分布的期望值μ决定了其位置，其标准差δ决定了分布的幅度。

当μ=0、δ=1时的正态分布是标准正态分布。10.1 Scipy中的科学计算工具3.统计分布【例10.5】正态分布应用。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('正态分布：$\mu$=%.1f,$\sigma^2$=%.1f'%(mu,sigma))plt.grid()plt.show()10.1 Scipy中的科学计算工具

4.积分egrate模块是SciPy中负责处理积分操作的模块。10.1 Scipy中的科学计算工具egrate模块中的常用函数函数说明quad一重积分dblquad二重积分tplquad三重积分nquadn倍多重积分fixed_quad高斯积分【例10.6】使用Segrate实现[0,1]区间上高斯函数的积分。importegratefromnumpyimportexpf=lambdax:exp(-x**2)i=egrate.quad(f,0,1)print(i)10.1 Scipy中的科学计算工具程序运行结果：(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中的科学计算工具程序运行结果：(0.375,1.2071790918144718e-14)5.插值10.1 Scipy中的科学计算工具插值是通过已知的离散点求未知数据的过程或方法。SciPy提供了erpolate模块进行插值操作：fromscipyimportinterpolate（1）一维插值

在给定有限数据点集合的情况下，通过构建一个函数来近似估计这些数据点之间的值。

interp1d()函数：一维数据的插值运算

该函数接收两个参数x点和y点，返回值是可调用函数

该函数可以用新的x调用并返回相应的y，y=f(x)。【例10.8】一维插值应用。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中的科学计算工具5.插值10.1 Scipy中的科学计算工具（2）样条插值在样条插值中，将点与由多项式定义的分段函数(称为样条)拟合。

InterpolatedUnivariateSpline()：样条插值函数

接受x和y作为参数并生成可调用的函数，该函数可通过新的x进行调用。【例10.9】样条插值应用。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中的科学计算工具

1.方程求解及求极值optimize.root()函数：求解非线性方程的根函数参数：

fun：表示方程的函数。

x0：根的初始猜测。10.2 Scipy中的优化10.2 Scipy中的优化1.方程求解及求极值optimize.minimize()函数：求函数最小值函数参数：

fun：要优化的函数

x0：初始猜测值

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

的根。fromscipy.optimizeimportrootfrommathimportcosdeffunc(x):returnx+cos(x)myroot=root(func,0)print(myroot.x)10.2 Scipy中的优化程序运行结果：[-0.73908513]

importnumpyasnpfromscipyimportoptimizeimportmatplotlib.pyplotaspltdeff(x):returnx**2+3*x+4x=np.arange(-10,10,0.1)plt.plot(x,f(x))#绘制函数曲线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中的优化程序运行结果：Optimizationterminatedsuccessfully.Currentfunctionvalue:1.750000Iterations:2Functionevaluations:6Gradientevaluations:3xmin:-1.5000000141634802ymin:1.7510.2 Scipy中的优化函数曲线及最小值2.数据拟合拟合是用一个连续函数（曲线）靠近给定的离散数据，使其与给定的数据相吻合。(1)多项式拟合用一个n阶多项式描述点(x,y)的关系：10.2 Scipy中的优化

多项式拟合的目的是找到一组系数a，使得拟合得到的曲线和真实数据点之间的距离最小。多项式拟合的系数可以使用Numpy模块的np.ployfit()函数得到。【例10.12】多项式拟合。importnumpyasnpimportmatplotlib.pyplotaspltx=np.linspace(-5,5,100)y=4*x+15#待拟合函数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中的优化2.数据拟合拟合结果如图10.7所示。10.2 Scipy中的优化2.数据拟合(2)曲线拟合曲线拟合可分为没有约束的曲线拟合和带有约束条件的曲线拟合。Scipy的optimize模块中提供了curve_fit()函数来拟合没有约束条件的曲线，least_squares()函数用来拟合带有约束条件的曲线。【例10.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中的优化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中的优化曲线拟合结果

10.3 Scipy中的稀疏矩阵处理在矩阵中，若数据为零的元素数目远远多于非零元素的数目，并且非零元素分布没有规律时，则称该矩阵为稀疏矩阵。根据存储方式不同，可以将稀疏矩阵分成以下几类：bsr_matrix：分块存储，基于行coo_matrix：坐标形式存储（COO）csr_matrix：基于行的压缩存储（CSR）csc_matrix：基于列的压缩存储（CSC）dia_matrix：对角线存储dok_matrix：基于键值对的存储lil_matrix：基于行的链表存储【例10.14】使用coo_matrix创建稀疏矩阵。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中的稀疏矩阵处理程序运行结果：

(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中的稀疏矩阵处理程序运行结果：array([[0,0,0,0],[8,0,1,7],[0,0,1,0],[0,0,0,0],[0,0,0,2]])【例10.16】使用data属性查看存储的数据（非零项）。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中的稀疏矩阵处理程序程序运行结果：[122]【例10.17】使用count_nonzero()方法计数非零元素的总数。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中的稀疏矩阵处理程序运行结果：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中的稀疏矩阵处理程序运行结果：

(1,0) 1(2,0) 2(2,2) 2【例10.19】使用sum_duplicates()方法消除重复项。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中的稀疏矩阵处理程序运行结果：

(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("相加结果：\n",(csrB+csrA).toarray())print("相减结果：\n",(csrB-csrA).toarray())print("相乘结果：\n",sparse.csc_matrix.multiply(csrA,csrB).toarray())10.3 Scipy中的稀疏矩阵处理程序运行结果：矩阵A：

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

[[6000][0012][0000][0010][0009]]相加结果：

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