Scipy学习


Scipy

address:https://scipy.github.io/devdocs/tutorial/index.html

Scipy内置模块

模块名 功能 参考文档
scipy.constants 数学常量 constants API
scipy.cluster 向量量化 cluster API
scipy.fftpack 快速傅里叶变换 fft API
scipy.integrate 积分 integrate API
scipy.interpolate 插值 interpolate API
scipy.io 数据输入输出 io API
scipy.linalg 线性代数 linalg API
scipy.ndimage N 维图像 ndimage API
scipy.odr 正交距离回归 odr API
scipy.optimize 优化算法 optimize API
scipy.signal 信号处理 signal API
scipy.sparse 稀疏矩阵 sparse API
scipy.spatial 空间数据结构和算法 spatial API
scipy.special 特殊数学函数 special API
scipy.stats 统计函数 stats.mstats API
scipy.misc 图像处理 misc API

Scipy常数模块

Scipy.constants模块内置多种常数,包括:

  • 公制单位
  • 字节为单位的二进制
  • 质量单位
  • 角度换算
  • 时间单位
  • 长度单位
  • 压强单位
  • 体积单位
  • 速度单位
  • 温度单位
  • 能量单位
  • 功率单位
  • 力学单位

Scipy cluster模块

这个模块分为vq模块和hierarchy模块.

vq模块

此模块实现k-means算法和向量量化.

K-means算法

  • 1.首先给定需要的类/簇个数n
  • 2.随机定n个簇心
  • 3.将每个点分在离它最近的簇心
  • 4.对于每个簇心,计算属于它的点集的平均点,该点为新的簇心
  • 5.在获得新的n个簇心之后,重复 3,4步骤,如果划分新的点集的时候,点集没有变则循环结束.

这个算法最后返回的是质心坐标矩阵,通过vq方法再把每个点分给对应的质心,vq方法的返回值是每个观察的簇是失真.

hierarchy模块

博客多元分析一节

Scipy fftpack模块

fftpack模块提供了快速傅里叶变化的方法.

一维离散傅里叶变换

长度为n的x序列,通过fft方法算出其在复数域上的映射a+bi.ifft方法是逆傅里叶变换.

离散余弦变换

离散余弦变换是实偶函数的傅里叶变换,不含虚数单位项.可以认为是0i.也即没有初项,只有振幅和频率两个参数.
Scipy提供了dct和idct方法求得

Scipy integrade模块

下表是integrate库中的常用功能

编号 功能说明
1 quad单积分
2 dblquad双重积分
3 tplquad三重积分
4 nquadn倍多重积分
5 fixed_quad高斯积分,阶数n
6 quadrature高斯正交到容差
7 romberg隆伯格积分
8 trapz梯形法则
9 cumtrapz梯形法则累计计算积分
10 simps辛普森的规则
11 romb隆伯格整合
12 polyint分析多项式积分(NumPy)
13 poly1d辅助函数polyint(NumPy)

quad 单积分

$\int_a^bf\left(x\right)dx$ == quad(f,a,b)

quad返回两个值,第一个是积分结果,第二个是积分的绝对误差

实例:

import scipy.integrate
from numpy import exp
f= lambda x:exp(-x**2)
i = scipy.integrate.quad(f, 0, 1)
print i

双重积分

对于积分
$$\int_a^bdx\int_{g(x)}^{h(x)}dyf(x,y)$$
其等价于dblquad(func,a,b,gfun,hfun)

其中func是f(x,y)这个复合函数 a,bx的上下限,gfunchfunc是y的上下限.

输出和单积分相同,一个是积分结果,一个是绝对误差,

interpolate 插值

什么是插值

插值就是在离散数据的基础上补插连续函数,使得连续函数通过这些离散点.这是离散函数逼近的重要方法.

Interp1d

scipy.interpolate中的interp1d类是一种基于固定数据点创建函数的便捷方法.

interp1d(datax,datay,kind)

其中的datax,datay分别为离散参数,kind是插值方式,比如线性(liner),最近(nearest),零,二次(cubic),立方等.

多维数据插值-griddata

griddata(points,value,point_grid,method)

参数解释:

  • points 是每个点的坐标
  • value 是点相应的函数结果值
  • point_grid 是你要想获得结果的点的坐标
  • method 是插值方法

这个方法返回的是点的函数结果

单变量样条

UnivariateSpline是基于固定数据点类创建函数的便捷方法.

scipy.interpolate.UnivariateSpline(x,y,w = None,bbox = [None,None],k = 3,s =无,ext = 0,check_finite = False)

参数解释:

  • ‘w’ - 指定样条拟合的权重。 必须是正数的。 如果没有(默认),权重都是相等的。
  • ‘s’ - 通过指定平滑条件指定结的数量。
  • ‘k’ - 平滑样条曲线的度数。 必须<= 5.默认值为k = 3,三次样条曲线。
  • Ext - 控制不在结节序列定义的区间内的元素的外推模式。
    • 如果ext = 0或’extrapolate’,则返回外推值。
    • 如果ext = 1或’0’,则返回0
    • 如果ext = 2或’raise’,则引发ValueError
    • 如果ext = 3’const’,则返回边界值。
  • check_finite - 是否检查输入数组是否仅包含有限数字。

Scipy I/O模块

MATLAB files

方法 描述
loadmat(file_name[, mdict, appendmat]) Load MATLAB file.
savemat(file_name, mdict[, appendmat, …]) Save a dictionary of names and arrays into a MATLAB-style .mat file.
whosmat(file_name[, appendmat]) List variables inside a MATLAB file.
  • loadmat

    squeeze_me参数为true时能将文件中1*1大小的矩阵滤除
    struct_as_record 参数为False时能将struct类型变成类似python的对象,能直接通过属性访问数据,而不是字典.

  • savemat 可以通过字典的形式存
    a_dict = {'field1': 0.5, 'field2': 'a string'}
    
    sio.savemat('saved_struct.mat', {'a_dict': a_dict})

    ID files

    方法 描述
    readsav(file_name[, idict, python_dict, …]) Read an IDL .sav file.

    Matrix Market files

    方法 描述
    mminfo(source) Return size and storage parameters from Matrix Market file-like ‘source’.
    mmread(source) Reads the contents of a Matrix Market file-like ‘source’ into a matrix.
    mmwrite(target, a[, comment, field, …]) Writes the sparse or dense array a to Matrix Market file-like target.

    Scipy.linalg 模块

    Scipy,linalg包含所有的numpy.linalg函数,另外还有一些高级函数.

    linalg.solve

    solve(a,b) a是系数矩阵,b是右手侧常数项.如
    #importing the scipy and numpy packages
    from scipy import linalg
    import numpy as np
    #Declaring the numpy arrays
    a = np.array([[3, 2, 0], [1, -1, 0], [0, 5, 1]])
    b = np.array([2, 4, -1])
    #Passing the values to the solve function
    x = linalg.solve(a, b)
    #printing the result array
    print x
    等于下列线性方程组
    $$3x + 2y + 0z = 2$$
    $$1x + -1y + 0z = 4$$
    $$0x + 5y + 1z = -1$$
    这个方法返回解数组

    行列式值

    Scipy.linalg.det函数计算矩阵的标量值

    特征值和特征向量

    $$Av = \lambda v$$
    A是矩阵,v是特征向量,$\lambda$是特征值.scipy.linalg.eig能返回函数特征值和特征向量.

l,v = eig(array) l是特征值,v是特征向量.

奇异值分解

理论解释

奇异值分解(SVD)可以认为是特征值问题对非平方矩阵的扩展.

对于方阵,如果求出了他的特征值和特征向量,那么方阵$A = W\Sigma W^T$

其中W为n个特征向量张成的矩阵,$\Sigma$是n个相应的特征值为主对角线的n*n维矩阵.

那么对于非平方矩阵,也有类似的特征分解表达式
$$A = U\Sigma V^T$$
A是一个m*n的矩阵,U是一个m*m的矩阵,$\Sigma$是一个m*n的矩阵,除了主对角线元素以外全为0,主对角线每个元素成为奇异值,V是一个n*n的矩阵.U和V都是酉矩阵,即满足$U^TU = I, V^TV = I$

那么$A^TA = V\Sigma^T U^TU\Sigma V^T =V\Sigma’ V^T = W\Sigma W^T$
所以由所有特征向量张成的矩阵就是矩阵V.

同样$AA^T = U\Sigma’’U^T$,那么其所有特征向量张成的矩阵就是矩阵U

对于奇异值矩阵我们发现$AV = U\Sigma$那么$Av_i=\sigma_iu_i$所以$\sigma_i=Av_i/u_i$

scipy.linalg里面的svd方法能求出相应的U,S,V矩阵.
U,S,V = svd(array)

Scipy.ndimage模块

ndimage模块主要用于图像处理,常见功能如下:

  • 输入/输出,显示图像
  • 基本操作 - 裁剪,翻转,旋转等
  • 图像过滤 - 去噪,锐化等
  • 图像分割 - 标记与不同对象相对应的像素
  • Classification
  • Feature extraction
  • Registration

    过滤器filter

    下面的函数都是对数组中的值进行过滤,输出元素是输入元素邻域元素的函数值.对于这个元素邻域我们成为过滤器内核,可以是矩形,也可以是任意图形,通过自己给的权重控制.对于过滤器内核,其内核中心一般是核形状的尺寸除以2,(长度为2的内核,中心就是1,既第一位元素.),但是也能通过设置参数移动内核中心.

对于边界元素,需要边界外的值来辅助计算,所以需要选择边界条件.

mode description example
nearest 使用边界值 [1,2]->[1,1,2,2]
wrap 复制阵列 [1,2,3]->[3,1,2,3,1]
reflect 反射在边界的数组 [1,2,3]->[1,1,2,3,3]
mirror 镜像边界处的数组 [1,2,3]->[2,1,2,3,2]
constant 使用常数填充(默认是0) [1,2,3]->[0,1,2,3,0]

同时还支持以下关键字

mode description
grid-constant 等价于constant
grid-mirror 等价于mirror
grid-wrap 等价于wrap
  • correlate是计算多维相关的函数,数组和给定内核相关

    scipy.ndimage.correlate(input, weights, output=None, mode='reflect', cval=0.0, origin=0)

    参数解释

    • input:输入的矩阵
    • weights:内核权重
    • output:要在其中放置输出的数组,或返回数组的 dtype。默认情况下,将创建与输入相同的 dtype 数组。
    • mode:对于边界元素的处理方式
    • cval:如果mode是constant,填充的常数,默认是0
    • origin: 内核中心偏移
  • correlate1d 是计算一维线性相关的函数

    scipy.ndimage.correlate1d(input, weights, axis=- 1, output=None, mode='reflect', cval=0.0, origin=0)[source]

    >from scipy.ndimage import correlate1d
    >correlate1d([2, 8, 0, 4, 1, 9, 9, 0], weights=[1, 3])
    >array([ 8, 26,  8, 12,  7, 28, 36,  9])

    插值函数

    Scipy.ndimage的插值函数是基于B-spline理论.其更详细的解释可以查看官方文档,本人因数理基础和时间原因,暂不能继续学习,未来或许会更新.

    形态学

    官方文档

    距离变换

    官方文档

    Scipy.odr正交距离回归

    odr表示Orthogonal Distance Regression,用于回归研究.在标准线性回归中,通常是通过计算Y的物产,但是有时候考虑X,Y和真实值的垂直距离误差更为明智.
    ODR方法需要传递data,线性回归模型.
    ODR(data,linear_model,beta0= [0,1])
    ## 参数解释
    linear_model = Model(linear_func)
    data = RealData(x,y)
    linear_func = lambda x,p:px+p
    x = array()
    y = array()

Scipy.Optimize

Optimize提供了几种优化算法,包括:

  • 无约束和约束最小化
  • 全局(强力)优化程序
  • 最小二乘最小化和曲线拟合
  • 标量单变量函数最小化器
  • 根查找
  • 使用各种算法的多变量方程系统求解器

    无约束和约束最小化多元标量函数

    minimize方法是多变量标量函数无约束或约束最小化算法的通用接口.不同算法需要不同的method. minimize至少需要提供三个参数,分别是函数,自变量数组,算法.

下面是rosenbrock函数的示例

import numpy as np
from scipy.optimize import minimize

def rosen(x):
    """The Rosenbrock function"""
    return sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0)

x0 = np.array([1.3, 0.7, 0.8, 1.9, 1.2])
res = minimize(rosen, x0, method='nelder-mead',
               options={'xatol': 1e-8, 'disp': True})

其他算法详解:https://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html

全局优化

有多个方法,分别是shgo,differential_evolution,basinhopping等,这些算法只需要提供函数和boundry就行.

最小二乘法优化

SciPy能够求解鲁棒约束的边界约束非线性最小二乘问题.给定残差f(x),损失函数rho(s),least_squares能找到代价函数F(x)的局部最小值

根查找

root方法能找到非线性方程的根.示例如下.

import numpy as np
 from scipy.optimize import root
def func(x):
...     return x + 2 * np.cos(x)
sol = root(func, 0.3"""猜测值""")
print(sol.x) //###array([-1.02986653])
print(sol.fun) //误差
###array([ -6.66133815e-16])

Scipy.sparse稀疏矩阵

稀疏矩阵有不同的表现形式,主要是

  • CSC -压缩稀疏列,按列压缩
  • CSR -压缩稀疏行,按行压缩

对于稀疏矩阵,有data方法查看其矩阵内部非零元素,count_nonzero()方法计算非零元素总数,对于稀疏矩阵运算和linalg模块类似,但是需要调用sparse的子模块,sparse.linalg模块内部函数,具体看文档.

scipy.spatial 空间数据结构算法

scipy.special 特殊数学函数

立方根函数

scipy.special.cbrt(x)获取x逐元素立方根。

指数函数

scipy.special.exp10(x)计算$10^x$。

scipy.special.exp2(x)计算$2^x$。

scipy.special.exp(x)计算$e^x$。

相对误差指数函数

scipy.special.exprel(x)生成相对误差指数,(exp(x) - 1)/ x。
当x接近零时,exp(x)接近1,因此exp(x)-1的数值计算可能遭受灾难性的精度损失。 所以需要xprel(x)以避免精度损失。

对数和指数函数

logsumexp(x) == log(sum(exp(a)))

朗伯W函数

lambertw(x)函数是$f(x) = x * e^x$的反函数,也就是说
$$lambertw(x) = w$$
$$x = w*e^w$$

排列组合

组合函数的语法是 - scipy.special.comb(N,k)等价于 $C_n^k$

排列函数的语法是 - scipy.special.perm(N,k)等价于 $A_n^k$

伽马功能

对于自然数n,伽马函数通常被称为广义阶乘,也即z * gamma(z)= gamma(z + 1)和gamma(n + 1)= n !,

scipy.stats 统计数学函数

文档


Author: Dovahkiin
Reprint policy: All articles in this blog are used except for special statements CC BY 4.0 reprint policy. If reproduced, please indicate source Dovahkiin !
  TOC