SciPy科学计算

常用SciPy子模块

模块	说明
cluster	向量量化
constants	数学常量
integrate	积分
interpolate	插值
optimize	优化
linalg	线性代数
stats	统计
io	数据的输入输出

SciPy中的科学计算工具

线性方程组求解

#scipy.linalg.solve(arrary1,arrary2)
import numpy as np
from scipy import linalg
a=np.array([[3,5,0],[1,-1,0],[0,3,1]])
b=np.array([5,3,-1])
x=linalg.solve(a,b)
print(x)

范数

范数(norm)是数学中的一种基本概念。在泛函分析中，它定义在赋范线性空间，并满足一定的条件，即1)非负；2)齐次性;3)三角不等式。它常常被用来度量某个向量空间（或矩阵）中的每个向量的长度或大小。
SciPy提供了函数linalg.norm()来计算向量或矩阵的范数，常用范数见下表，

向量范数	含义	矩阵范数	含义
0范数	向量中非零元素个数	1范数	矩阵中每列元素和的最大值
1范数	向量元素绝对值之和	-1范数	矩阵中每列元素和的最小值
2范数	欧几里得范数	'fro’范数	矩阵中所有元素平方和的平方根
∞范数	所有向量元素绝对值中的最大值	∞范数	矩阵中每行元素和的最大值
-∞范数	所有向量元素绝对值中的最小值	-∞范数	矩阵中每行元素和的最小值
p范数	向量元素绝对值的p次方和的1/p次幂	2范数	矩阵的最大奇异值

from scipy import linalg
import numpy as np
M=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'))

统计分布

SciPy中的stats模块包含了各种连续分布和离散分布的模型。在统计分布中，用概率质量函数（Probability Mass Function，PMF）和累积分布函数（Cumulative Distribution Function，CDF）表示离散变量的分布，用概率密度函数（Cumulative Distribution Function，PDF）和累积分布函数（Probability Distribution Function）表示连续变量的分布。PMF 定义了随机变量的结果小于或等于值 X 的概率，PDF 和 PMF 相同，但用于连续值，CDF 表示随机变量 X 的所有可能值 x 的概率，CDF 既可用于离散分布，也用于连续分布。

1. 伯努利分布
   伯努利分布是一种离散分布，又称为两点分布或 0-1 分布，它的随机变量只取 0 或者 1。进行一次事件试验，该事件发生的概率为 p，不发生的概率为 1-p。这是一个最简单的分布，任何一种只有两种结果的随机现象都服从于 0-1 分布，如抛硬币观察正反面、打枪上靶统计等。
   例子：打枪上靶实验。

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
import matplotlib as mpl
#设置图例字体
mpl.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')
#vline用于绘制竖直线(x坐标值,y坐标最小值,y坐标最大值)
plt.vlines(x,(0,0),pList)
plt.xlabel('随机变量：打枪上靶一次')
plt.ylabel('概率')
plt.title('伯努利分布:p=%.2f'%p)
plt.show()

2. 二项分布
   二项分布又称作n重伯努利分布。如果随机变量序列Xn(n=1,2,…)中的随机变量均服从参数为p的伯努利分布，那么随机变量序列Xi就形成了参数为p的n重伯努利分布。例如，假定重复抛掷一枚均匀硬币n次，如果在第i次抛掷中出现正面，令Xi=1；如果出现反面，则令Xi=0。那么，随机变量Xn(n=1,2,…)就形成了参数为1/2的n重伯努利试验。

例子：打枪上靶实验，打靶10次，命中8环的概率分布。

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
import matplotlib as mpl
#设置图例字体
mpl.rcParams['font.sans-serif'] = ['SimHei']
plt.grid()
n=10 #打枪次数
p=0.75 #8环命中的概率
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=%.2f"%(n,p))
plt.show()

3. 正态分布
   正态分布也称为高斯分布，是统计学中最主要的连续概率分布。
   若随机变量 $X$ 服从一个数学期望为 $\mu$、方差为 ${\delta }^2$ 的正态分布，记为 $N(\mu ,{\delta }^2)$。其概率密度函数为正态分布的期望值 $\mu$ 决定了其位置，其标准差 $\sigma$ 决定了分布的幅度。当 $\mu =0$ ，$\delta =1$ 时的正态分布是标准正态分布。

例子，正态分布应用。

import numpy as np
from scipy import stats
import matplotlib as mpl
import matplotlib.pyplot as plt
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus']=False
mu, sigma = 0, 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('标准正态分布')
plt.grid()
plt.show()

积分

scipy.integrate模块是SciPy中负责处理积分操作的模块，该模块中常用函数见下表

函数	说明
quad()	一重积分
dblquad()	二重积分
tplquad()	三重积分
nquad()	n倍多重积分
fixed_quad()	高斯积分
例子：使用scipy.integrate实现[0，1]区间上高斯函数的积分${\int }_0^1{e}^{-x^2}dx$。

import scipy.integrate
from numpy import exp
f=lambda x:exp(-x**2)
i=scipy.integrate.quad(f,0,1)
#使用quad()函数计算积分，该函数返回两个值，其中第一个是积值，第二个是积值绝对误差的估计值。
print(i)

例子：使用scipy.integrate计算二重积分${\int }_0^{1}dy{\int }_0^{\sqrt{1-y^2}}8xydx$。

import scipy.integrate
from math import sqrt
f=lambda x,y:8*x*y
g=lambda x:0
h=lambda y:sqrt(1-2*y**2)
i=scipy.integrate.dblquad(f,0,0.5,g,h)
print(i)

插值

插值是通过已知的离散点求未知数据的过程或方法。Scipy提供了scipy.interpolate模块进行插值操作，使用该模块时应先导入，可采用以下方法导入：

from scipy import interpolate

1. 一维插值

一维插值是一种在给定有限数据点集合的情况下，通过构建一个函数来近似估计这些数据点之间的值。一维数据的插值运算可以通过interp1d()函数完成，该函数接收两个参数x点和y点，返回值是可调用函数，该函数可以用新的x调用并返回相应的y，y=f(x)。

例子

from scipy.interpolate import interp1d
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(200)
x=np.random.randn(10)
y=np.random.randn(10)
f=interp1d(x,y,kind='cubic')
x_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()

2. 样条插值
   在一维插值中，将点拟合为一条曲线，在样条插值中，将点与由多项式定义的分段函数（称为样条）拟合。样条插值使用InterpolatedUnivariateSpline()函数，InterpolatedUnivariateSpline()函数接收x和y作为参数并生成可调用的函数，该函数可通过新的x进行调用。

例子

from scipy.interpolate import InterpolatedUnivariateSpline
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(10000)
x=np.linspace(-5,5,20)
y=np.exp(-x**2)+0.01*np.random.randn(20)+1
f=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()

SciPy中的优化

方程求解及求极值

SciPy 中的 optimize.root() 函数可以求解非线性方程的根，该函数接收两个参数，分别是：fun：表示方程的函数。x0：根的初始猜测。
该函数返回一个解决方案表示为 OptimizeResult 对象，其重要属性是：x 为解决方案的数组，success 为指示算法是否成功退出的布尔标志和说明终止原因的信息。
例子：使用 root() 函数求解方程 x + cos(x) = 0 的根。

from scipy.optimize import root
from math import cos
def func(x):
    return x+cos(x)
myroot=root(func,0)
print(myroot.x)

Scipy 中的 optimize.minimize() 函数用来求函数最小值，该函数接收三个参数，分别是：fun：要优化的函数。x0：初始猜测值。method：“BFGS”，“Newton-CG”，“L-BFGS-B”，“TNC”，“COBYLA”，“SLSQP”。
例子：使用 BFGS() 函数抛物线 f(x)=x^2+3x+4 的极小值。

import numpy as np
from scipy import optimize
import matplotlib.pyplot as plt
def f(x):
    return x**2 + 3*x + 4
x = np.arange(-10, 10, 0.1) # 绘制函数曲线
plt.plot(x, f(x))
# 第一个参数是函数名，第二个参数是梯度下降的起点，返回值是函数最小值的 x 值
xout = optimize.fmin_bfgs(f, 0)
xmin = xout[0] # x 值，即函数最小值
ymin = f(xmin) # y 值，即函数最小值
print('xmin:', xmin)
print('ymin:', ymin)
plt.scatter(xmin, ymin, s=15, c='r') # 绘制最小值的点
plt.show()

数据拟合

数据拟合是用一个连续函数（曲线）靠近给定的离散数据，使其与给定的数据相吻合。

1. 多项式拟合
   多项式拟合是一个n阶多项式描述(x,y)的关系：
   $y=a_n{x}^n+a_{n-1}{x}^{n-1}+\cdots +a_1{x}^1+a_0{x}^0={\sum }_{i=0}^n{a}_i{x}^i$

多项式拟合的目的找到一组系数a，使得拟合得到的曲线和真实数据点之间的距离最小。
多项式系数可以使用NumPy模块的np.polyfit()函数得到。

例子：多项式拟合。

import numpy as np
import matplotlib.pyplot as plt
x=np.linspace(-5,5,100)
y=4*x+15 # 待拟合函数
y_noise=y+np.random.randn(100)*2 # 添加随机噪声
coeff=np.polyfit(x,y_noise,1)
plt.plot(x,y_noise,'*r',x,coeff[0]*x+coeff[1])
plt.show()

2. 曲线拟合

曲线拟合可分为没有约束条件的曲线拟合和带有约束条件的曲线拟合。SciPy的optimize模块中提供了curve_fit()函数来拟合没有约束条件的曲线，least——squares()函数来拟合带有约束条件的曲线。

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
np.random.seed(10)
def fun(x):
    #用来模拟测量误差
    measurement_error = np.random.rand(len(x)) #产生随机数
    y = 5 * x ** 2 + measurement_error
    return y
def objective_fun(x,a):
    #定义二次函数作为目标函数
    return a * x**2
x = np.arange(0,10,0.2)
y = fun(x)
pout,m = curve_fit(objective_fun,x,y) #获取参数
a = pout
#绘制实际数据和预测的函数
plt.scatter(x,y,label='data')
plt.plot(x,objective_fun(x,a),"--",color="r",label="objective quadratic function")
plt.legend()
plt.show()

SciPy中的稀疏矩阵处理

在矩阵中，若数据为零的元素数目远远多于非零元素的数目，并且非零元素分布没有规律时，则称该矩阵为稀疏矩阵；与之相反的称为稠密矩阵。在科学与工程领域中求解线性模型时经常出现大型的稀疏矩阵。

稀疏矩阵的存储

对于稀疏矩阵，如果按照正常方式存储所有元素则这些矩阵所占空间将十分巨大，尤其是当矩阵的维度特别多的时候，稀疏矩阵会造成巨大的空间浪费。为了降低存储空间，稀疏矩阵通常只保存非零值及其对应的位置。

在Python中，根据存储方式不同，可以将稀疏矩阵分成以下几类：

bsr_matrix (Block Sparse Row matrix):基于行的分块存储。
coo_matrix (A sparse matrix in COOrdinate format):坐标形式存储(COO)。
csr_matrix (Compressed Sparse Row matrix):基于行的压缩存储(CSR)。
csc_matrix (Compressed Sparse Column matrix):基于列的压缩存储(CSC)。
dia_matrix (Sparse matrix with DIAGonal storage):基于对角线存储。
dok_matrix (Dictionary Of Keys based sparse matrix):基于键值对的存储。
lil_matrix (Row-based linked list sparse matrix):基于行的链表存储。

在存储方式中，COO方式在构建矩阵时比较高效，CSC 和 CSR 方式在乘法计算中较高效。
(1) Coordinate Matrix 对角存储矩阵(COO) coo_matrix是最简单的存储方式，采用三个数组 row、col 和 data 保存非零元素的行下标、列下标与值，这三个数组的长度相同。一般来说，coo_matrix 主要用来创建矩阵，因为 coo_matrix 无法对矩阵的元素进行增删改等操作，一旦创建之后，除了将其转换成其他格式的矩阵，几乎无法对其做任何操作和矩阵运算。

（2）csr_matrix 基于行的压缩存储（CSR
csr_matrix 是按照行对矩阵进行压缩的，通过 indices, indptr, data 三个数组来确定矩阵。data 表示矩阵中的非零数据，第 i 行非零元素的列索引为 indices[indptr[i]:indptr[i+1]]，可以将 indptr 理解成利用其自身索引来指向第 i 行元素的列索引，根据 [indptr[i]:indptr[i+1]]，可以得到该行中的非零元素个数。
若 index[i] = 5 且 index[i+1] = 5，则第 i 行没有非零元素，若 index[j] = 6 且 index[j+1] = 8，则 j 行的非零元素的列索引为 indices[6:8]，得到了行索引、列索引，相应的数据存放在：data[index[i]:index[i+1]]。

（3）csc_matrix 基于列的压缩存储（CSC）
csc_matrix 是按列对矩阵进行压缩的，通过 indices, indptr, data 三个数组来确定矩阵，与 CSR 类似，data 表示矩阵中的非零数据，第 i 列非零元素的行索引为 indices[indptr[i]:indptr[i+1]]，indptr 可以理解为利用其自身索引来指向第 i 列元素的列索引，根据 [indptr[i]:indptr[i+1]]，可以得到该列中的非零元素个数。
若 index[i] = 5 且 index[i+1] = 5，则第 i 列没有非零元素，若 index[j] = 6 且 index[j+1] = 8，则第 j 列的非零元素的行索引为 indices[6:8]，得到了列索引和行索引，相应的数据存放在：data[index[i]:index[i+1]]。

例子：使用 coo_matrix 创建稀疏矩阵。

import scipy.sparse as sp

values = [3, 2, 1, 0]
rows = [0, 1, 2, 3]
cols = [1, 3, 2, 0]
A = sp.coo_matrix((values, (rows, cols)), shape=[4, 4]) # 创建稀疏矩阵
print(A)
A.todense() # 将稀疏矩阵对象转为NumPy 的 array

例子：存储稀疏矩阵。

import numpy as np
import scipy as sp
indptr = 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 = sp.sparse.csr_matrix((data, indices, indptr))
csr.toarray()
print(csr)

例子：使用 data 属性性查看存储的数据（非零项）。

import numpy as np
from scipy.sparse import csr_matrix
arr = np.array([[0, 0, 0, 0], [1, 0, 0, 0], [2, 0, 2, 0], [0, 0, 0, 0]])
print(csr_matrix(arr).data)

例子：使用 count_nonzero() 方法统计非零元素的总数。

import numpy as np
from scipy.sparse import csr_matrix
arr = np.array([[0, 0, 0, 0], [1, 0, 0, 0], [2, 0, 2, 0], [0, 0, 0, 0]])
print(csr_matrix(arr).count_nonzero())#统计非零元素的个数

例子：使用 eliminate_zeros() 方法从矩阵中删除零项。

import numpy as np
from scipy.sparse import csr_matrix
arr = 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)

例子：使用 sum_duplicates() 方法消除重复项。

import numpy as np
from scipy.sparse import csr_matrix
arr = np.array([[0, 0, 0, 0], [1, 0, 0, 0], [2, 0, 2, 0], [0, 0, 0, 0]])
A = csr_matrix(arr)
A.sum_duplicates()
print(A)

稀疏矩阵的运算

Scipy 中的 sparse 提供了许多稀疏矩阵的运算，包括矩阵加法、减法、乘法以及矩阵求逆等。 在例10.20稀疏矩阵加法、减和乘法运算。
例子：稀疏矩阵加法、减法和乘法运算.

import numpy as np
from scipy import sparse
indptr=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())
indptrl=np.array([0,1,3,3,4,5])
indices1=np.array([0,2,3,2,3])
datal=np.array([6,1,2,1,9])
csrB=sparse.csr_matrix((datal,indices1,indptrl))
print("矩阵B:\n",csrB.toarray())
print("相加结果:\n",(csrB+csrA).toarray())# 加法运算
print("相减结果:\n",(csrB-csrA).toarray())# 减法运算
print("相乘结果:\n",sparse.csc_matrix.multiply(csrA,csrB).toarray())#乘法运算