SciPy简单应用

c = StringIO(u"1,0,2\n3,0,4")
"""
delimiter指定文件中数据之间的分隔符
usecols 指定读取文件中的哪些列，在下面的函数中，usecols=(0,2)表示只读取文件的第1、3列
unpack:读出来的数据是否需要解包,注意接受函数的变量个数要与读出来的行数或列数吻合
"""
x, y = np.loadtxt(c, delimiter=',', usecols=(0, 2), unpack=True)
x

array([1., 3.])

y

array([2., 4.])

特殊函数：scipy.special

常用的一些函数如下：

• 贝塞尔函数，比如scipy.special.jn() (第n个整型顺序的贝塞尔函数)
• 椭圆函数 (scipy.special.ellipj() Jacobian椭圆函数, …)
• Gamma 函数: scipy.special.gamma(), 也要注意 scipy.special.gammaln() 将给出更高准确数值的 Gamma的log。
• Erf, 高斯曲线的面积：scipy.special.erf(),也叫误差函数

线性代数操作：scipy.linalg

• scipy.linalg.det() 函数计算方阵的行列式：

from scipy import linalg
arr = np.array([[1, 2],[3, 4]])
linalg.det(arr)

-2.0

arr = np.array([[3, 2],[6, 4]])
linalg.det(arr) # 科学计数法 6.661338147750939e-16 =6.661338147750939 x 10^-16 约等于0

6.661338147750939e-16

• scipy.linalg.inv() 函数计算逆方阵

arr = np.array([[1, 2],[3, 4]])
iarr = linalg.inv(arr)
iarr

array([[-2. ,  1. ],[ 1.5, -0.5]])

np.allclose(np.dot(arr, iarr), np.eye(2)) # numpy的allclose方法，比较两个array是不是每一元素都相等，默认在1e-05的误差范围内

True

计算行列式为0的矩阵将抛出LinAlgError :

arr = np.array([[3, 2],[6, 4]])
linalg.inv(arr)

---------------------------------------------------------------------------LinAlgError                               Traceback (most recent call last) in 1 arr = np.array([[3, 2],2                  [6, 4]])
----> 3 linalg.inv(arr)~\Anaconda3\envs\tf2\lib\site-packages\scipy\linalg\basic.py in inv(a, overwrite_a, check_finite)972         inv_a, info = getri(lu, piv, lwork=lwork, overwrite_lu=1)973     if info > 0:
--> 974         raise LinAlgError("singular matrix")975     if info < 0:976         raise ValueError('illegal value in %d-th argument of internal 'LinAlgError: singular matrix

• 奇异值分解（SVD）
  A = U Σ V T A=U\varSigma V^T A=UΣVT

arr = np.arange(9).reshape((3, 3)) + np.diag([1, 0, 1])
print(arr)
uarr, spec, vharr = linalg.svd(arr)

[[1 1 2][3 4 5][6 7 9]]

结果的数组频谱(矩阵的特征值$\lambda $组成的对角阵)是：

spec

array([14.88982544,  0.45294236,  0.29654967])

U矩阵为：

uarr

array([[-0.1617463 , -0.98659196,  0.02178164],[-0.47456365,  0.09711667,  0.87484724],[-0.86523261,  0.13116653, -0.48390895]])

V矩阵为：

vharr

array([[-0.45513179, -0.54511245, -0.70406496],[ 0.20258033,  0.70658087, -0.67801525],[-0.86707339,  0.45121601,  0.21115836]])

根据公式，我们可以根据UV矩阵和特征值矩阵还原原始矩阵

sarr = np.diag(spec) # 将spec矩阵化为对角阵
svt_mat = uarr.dot(sarr).dot(vharr)
print(svt_mat)

[[1. 1. 2.][3. 4. 5.][6. 7. 9.]]

• 快速傅立叶变换：scipy.fftpack

假设一个（有噪音）的信号输入是这样：

time_step = 0.02
period = 5.
time_vec = np.arange(0, 20, time_step)# sin(wt+n)
sig = np.sin(2 * np.pi / period * time_vec) + \0.5 * np.random.randn(time_vec.size)

假设信号来自真实的函数，因此傅立叶变换将是对称的。scipy.fftpack.fftfreq() 函数将生成样本序列，而将计算快速傅立叶变换：

from scipy import fftpack
sample_freq = fftpack.fftfreq(sig.size, d=time_step) # 通过DFFT生成样本序列
sig_fft = fftpack.fft(sig)

因为生成的幂是对称的，寻找频率只需要使用频谱为正的部分：

pidxs = np.where(sample_freq > 0)
freqs = sample_freq[pidxs]
power = np.abs(sig_fft)[pidxs]

寻找信号频率：

freq = freqs[power.argmax()] # 频谱中功率的最高点就是信号的原始频率
np.allclose(freq, 1./period)  # 检查是否找到了正确的频率

True

现在高频噪音将从傅立叶转换过的信号移除：

sig_fft[np.abs(sample_freq) > freq] = 0

生成的过滤过的信号可以用scipy.fftpack.ifft()函数做傅里叶逆变换：

main_sig = fftpack.ifft(sig_fft)

查看结果：

plt.figure()
plt.plot(time_vec, sig )
plt.plot(time_vec, np.real(main_sig), linewidth=3 )
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.show()

优化及拟合：scipy.optimize

• 寻找标量函数的最小值

from scipy import optimize

让我们定义下面的函数：

def f(x):return x**2 + 10*np.sin(x)x = np.arange(-10, 10, 0.1)
plt.plot(x, f(x))
plt.show()

这个函数在-1.3附近有一个全局最小并且在3.8有一个局部最小。
找到这个函数的最小值的常用有效方式是从给定的初始点开始进行一个梯度下降。BFGS算法是这样做的较好方式：

optimize.fmin_bfgs(f, 0)

Optimization terminated successfully.Current function value: -7.945823Iterations: 5Function evaluations: 18Gradient evaluations: 6array([-1.30644012])

这个方法的一个可能问题是，如果这个函数有一些局部最低点，算法可能找到这些局部最低点而不是全局最低点，这取决于初始点：

optimize.fmin_bfgs(f, 3, disp=0) # 将起点设置为3，无法找到全局最低点

array([3.83746709])

如果我们不知道全局最低点，并且使用其临近点来作为初始点，那么我们需要付出昂贵的代价来获得全局最优。要找到全局最优点，最简单的算法是暴力算法，算法中会评估给定网格内的每一个点：

grid = (-10, 10, 0.1)
xmin_global = optimize.brute(f, (grid,))
xmin_global

array([-1.30641113])

• 寻找标量函数的根scipy.optimize.fsolve()：

要寻找上面函数f的根，比如f(x)=0的一个点

root = optimize.fsolve(f, 1)  # 我们的最初猜想是1
root

array([0.])

注意只找到一个根。检查f的图发现在-2.5左右还有应该有第二个根。通过调整我们最初的猜想，我们可以发现正确的值：

root2 = optimize.fsolve(f, -2.5)
root2

array([-2.47948183])

• 曲线拟合

假设我们有来自f的样例数据，带有一些噪音：
y=h(t) * x(t)+n(t)

xdata = np.linspace(-10, 10, num=20)
ydata = f(xdata) + np.random.randn(xdata.size) # 第二项为噪音

现在，如果我们知道这些sample数据来自的函数（这个案例中是 x 2 + s i n ( x ) x^2 + sin(x) x2+sin(x)）的函数形式，而不知道每个数据项的系数，那么我们可以用最小二乘曲线拟合在找到这些系数。首先，我们需要定义函数来拟合：

def f2(x, a, b): # a,b为拟合参数，即x^2 和sin(x) 的系数return a*x**2 + b*np.sin(x)

然后我们可以使用scipy.optimize.curve_fit()来找到a和b：

guess = [2, 2] # 猜测的参数，合理设置该参数可以大幅度降低计算量
params, params_covariance = optimize.curve_fit(f2, xdata, ydata, guess)
params

array([ 0.99987461, 10.45777129])

plt.plot(xdata,ydata,linewidth=3)
a,b = params
plt.plot(xdata,f2(xdata,a,b))

[]

统计和随机数：scipy.stats

常见分布

可能用到的分布对照表

通用函数

stats连续型随机变量的公共方法：

• 离散分布的简单方法大多数与连续分布很类似，但是pdf被更换为密度函数pmf。

import numpy as np
from scipy import stats

a = np.random.normal(size=1000) # 生成正态分布
bins = np.arange(-4, 5) # -4~5的采样点
# a
bins

array([-4, -3, -2, -1,  0,  1,  2,  3,  4])

histogram = np.histogram(a, bins=bins, normed=True)[0]
bins = 0.5*(bins[1:] + bins[:-1])
bins

C:\Users\ASUS\Anaconda3\envs\tf2\lib\site-packages\ipykernel_launcher.py:1: VisibleDeprecationWarning: Passing `normed=True` on non-uniform bins has always been broken, and computes neither the probability density function nor the probability mass function. The result is only correct if the bins are uniform, when density=True will produce the same result anyway. The argument will be removed in a future version of numpy."""Entry point for launching an IPython kernel.array([-3.5, -2.5, -1.5, -0.5,  0.5,  1.5,  2.5,  3.5])

from scipy import stats
import pylab as pl
b = stats.norm.pdf(bins)  # 相当于np.random.normal，生成正态分布
pl.plot(bins, histogram)
pl.plot(bins, b)
pl.show()

• 生成服从指定分布的随机数

norm.rvs通过loc和scale参数可以指定随机变量的偏移和缩放参数，这里对应的是正态分布的期望和标准差。size得到随机数数组的形状参数。(也可以使用np.random.normal(loc=0.0, scale=1.0, size=None))

normpdf = stats.norm.rvs(loc=0,scale=1.0,size=10)
normpdf

array([-0.35670277, -0.35593883,  0.55816991, -1.02267389, -0.26461898,1.46413836, -1.11843986, -1.67888616,  0.39230616, -0.32538297])

np.random.normal(loc=0.0, scale=1.0, size=10)

array([-0.17223573,  0.34488558,  1.30292964, -0.29186785, -0.35257427,-1.99991253, -0.27301159,  1.52191222,  0.51331207,  0.77448018])

stats.norm.rvs(loc = 3,scale = 10,size=(2,2))

array([[-8.97012362,  0.82079641],[ 6.5068092 ,  8.60501163]])

• 求概率密度函数指定点的函数值

# stats.norm.pdf正态分布概率密度函数。
stats.norm.pdf(0,loc = 0,scale = 1) # 求原点的函数值

0.3989422804014327

stats.norm.pdf(np.arange(3),loc = 0,scale = 1)# 求0,1,2三个点的服从N~(0,1)标准正态分布的函数值

array([0.39894228, 0.24197072, 0.05399097])

本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！