4

SciPy库¶

1. SciPy是一款方便、易于使用、专为科学和工程设计的Python工具包。2. 它包括统计,优化,整合,线性代数,傅里叶变换,信号和图像处理,常微分方程求解器等等模块。3. SciPy库为很多问题提供了便捷且高效的数值路径解法。

---

要说明的是，在导入SciPy模块时，一般针对想实现的功能导入相应的子模块即可，不需要导入完整的SciPy模块。
在本章中通过例子的形式来讲解某几个子模块下(表格中加粗显示)，需要说明的是这里注重的是模块的应用，具体的算法或者原理只做简单陈述

SciPy子模块构成¶

子模块	功能	子模块	功能
cluster	聚类算法	**optimize	优化和求根**
constants	物理数学常数	signal	信号处理
fftpack	快速傅里叶变换	sparse	稀疏矩阵
integrate	积分和常微分方程求解	spatial	空间数据结构和算法
interpolate	插值处理	special	特殊方程
io	输入输出	stats	统计分布和函数
linalg	线性代数	weave	C/C++积分
odr	正交距离回归

​x

import numpy as np

一、积分和微分方程求解——integrate¶

利用scipy模块下的integrate子模块可以进行积分求解，和微分方程的求解。

1.1 数值积分—quad¶

quad函数求积分，会输出积分值，和最大误差。注意quad函数一次只能求解一重积分。$$\int_{a}^{b}dF(x) = ?$$$$\int_{a}^{b}dF(x) = lim\sum_{i=1}^{n}f(x_i) \Delta x_i$$

1

from scipy.integrate import quad #积分求解

1

x = np.linspace(0, 10, 50)

upper bound on error: 7.2e-11

1.2 常微分方程求解—odeint¶

$$\frac{dy}{dx}=sin(x)$$

1

def fx(y,x):

1

from scipy.integrate import odeint

[[ 0. ] [ 0.00201328] [ 0.00804512] [ 0.01807123] [ 0.03205123] [ 0.04992881] [ 0.071632 ] [ 0.0970734 ] [ 0.12615056] [ 0.1587464 ]]

1

plt.figure(figsize=(12,4))

二、最优化求解——optimize¶

2.1 曲线拟合¶

曲线拟合实际上是一个最优化某目标函数的问题

1

from scipy.stats import norm #norm用来产生正态随机数 在第三部分会详细讲解

1

def function(x, a , b, f, phi):

人为构造数据集，真实函数如下图红线所示，蓝色点表示的是样本点

1

x = np.linspace(0, 2 * np.pi, 50)

假设我们已经知道了函数的形式为$y = ae^{-bsin(fx+\phi)}$，只需要估计函数中的参数
基于样本数据，利用最小二乘的方法来估计函数中的参数

1

from scipy.optimize import leastsq

1

def f_err(p, y, x):

1

plt.plot(x,y, label='真实函数')

2.2 资产组合的有效边界求解¶

例如我们有四种资产，利用这四种资产构建投资组合，根据马克维茨的现代投资组合理论可以求解有效边界。求解有效边界的过程实际上就是一个求解最优化问题的过程。

1

r = np.array([-0.15, 0.12, 0.31, 0.04])

为了更好突显有效边界，这里随机构造4000个投资组合

1

port_r, port_var = [], []

利用马克维茨的MPT——在组合目标收益率情况下，最小化组合方差，求解有效边界。具体而言就是在不同的期望回报率约束下，求解下面的数学优化问题

$$Min 　投资组合方差=\sum_{i=1}^{n}\sum_{j=1}^{n}w_iw_j\sigma_{ij}$$
$$s.t. \sum_{i=1}^{n}w_iE(r_i) = 期望回报率, \sum_{i=1}^{n}w_i = 1$$

1

import scipy.optimize as sco #求解约束最优化问题

1

def statistics(weights):

可视化有效边界，可看到随机构造的投资组合都越不过有效边界

1

risk_free = 0.04

三、概率统计¶

3.1 描述性统计¶

1

import scipy.stats as stats

1

x_sample = np.random.randint(0,100,20)

[21 31 0 30 57 29 98 39 5 6 73 81 81 90 29 79 46 46 23 10]

1

print('样本均值：',x_sample.mean())

样本均值： 43.7样本中位数： 35.0样本最小值： 0样本最大值： 98样本总和： 874样本标准差： 29.8983277124样本众数： ModeResult(mode=array([29]), count=array([2]))样本2阶矩： 893.91样本偏度： 0.31757144674066823样本峰度： -1.1675337430657668

3.2 连续概率分布¶

scipy的统计模块下，包含了很多种连续概率分布函数。对于这些连续概率分布函数，通常会用四种操作：

• 分布简称.cdf——累计分布函数
• 分布简称.pdf——概率密度函数
• 分布简称.rvs——指定参数的随机变量
• 分布简称.fit——对于某个样本拟合该分布，对各个参数的MLE值

3.2.1 正态分布¶

1

from scipy.stats import norm

1

x = np.linspace(-3,3,50)

1

x_norm = norm.rvs(0,4,100)

均值0标准差4的正态随机数10个： [ -2.8629 -3.1896 3.7239 8.5456 -0.4123 -10.5606 -5.7586 2.112 3.258 -0.7308 5.5889 -2.0959 -3.1837 5.4953 -0.8358 1.1265 -1.1252 -1.2015 6.2171 -2.736 4.4792 -2.1679 4.4727 1.5076 -3.0379 -0.4054 7.2151 -5.6364 1.4001 0.897 -1.5876 1.8073 4.0069 -2.3222 -0. -0.64 0.3895 0.8823 -3.0253 -2.1211 3.7191 -3.6372 5.0046 0.5788 5.693 1.2093 6.468 2.5056 -2.9798 7.2322 -0.4794 -2.9561 -4.5213 -3.488 7.2058 4.8799 7.6271 -5.6004 1.0138 5.0218 3.0625 -7.9592 -0.1816 2.7152 2.0945 -0.4867 -3.7178 -1.0189 -4.1217 -4.7124 0.3075 -9.1469 8.3258 0.7299 0.3687 0.0684 -1.0141 3.2905 2.1668 -4.2256 -8.2148 -4.242 2.9949 -1.5833 -4.5376 5.1325 4.7154 -4.5276 2.9836 3.6237 -0.1887 -7.3434 -0.5845 -0.4368 7.502 0.2786 -2.1179 3.124 -5.3939 0.7396]

1

plt.hist(x_norm, normed=True, bins=20,color='blue',label='正态样本分布图')

1

x_mean, x_std = norm.fit(x_norm)

基于样本x_norm的总体均值的MLE： 0.20453974482基于样本x_norm的总体标准差的MLE： 4.11418088811

通过之前的integrate模块的积分功能可以帮助求解概率——概率密度函数的积分

1

x1 = np.linspace(-1,2,100)

1

from scipy.integrate import trapz

1

p = trapz(norm.pdf(x1), x1)

81.86% of the values lie between -1 and 2

3.2.2 其他分布¶

1

from scipy.stats import lognorm, chi2, f, t #导入对数正态分布 卡方分布 F分布 t分布

对数正态分布

1

x = np.linspace(0.01, 3, 100)

卡方分布

1

x = np.linspace(0.01, 30, 100)

F分布

1

x = np.linspace(0.01, 10, 100)

t分布

1

x = np.linspace(-3, 3, 100)

3.3 离散概率分布¶

对于离散概率分布函数，通常会用三种操作：

• 分布简称.cdf——累计分布函数
• 分布简称.pdf——概率质量函数
• 分布简称.rvs——指定参数的随机变量

1

from scipy.stats import randint, binom, poisson #离散均匀分布，二项分布，泊松分布

离散均匀分布

1

x = np.arange(0, 11, 1)

1

x = np.arange(0, 11, 1)

1

x_dis_uniform = randint(0,10).rvs(10)

产生10个0-10离散均匀随机数： [0 6 4 9 1 7 8 0 6 4]

二项分布

1

x = np.arange(0, 50, 1)

泊松分布

1

x = np.arange(0, 20, 1)

3.4 假设检验¶

3.4.1 正态性检验¶

检验一组样本是否可认为服从正态分布。例如我们想检验下面的x_for_test这个样本是否可认为来自正态分布的总体。

1

x_for_test = np.linspace(-20,20,20)

KS检验，从峰度和偏度出发，进行的检验。该检验不光能检验某数据集服从正态，也可以检验该数据集是否服从其他的分布。
用法：kstest(待检验数据集, 何种类型分布的检验, args=(), N=20, alternative=’two_sided’(设置双尾，单尾检验等), mode=’approx’, **kwds)
具体的可以参考：https://docs.scipy.org/doc/scipy-0.7.x/reference/generated/scipy.stats.kstest.html

1

from scipy.stats import kstest

经过KS检验发现p值很小，则拒绝原假设，即有充足理由认为该样本不是来自于正态总体的

Shapiro检验，专门用来做正态性检验的模块。需注意的是，在样本数很大的情况下，该检验不适合做正态性检验，检验结果可能不准确。
具体可以参考：https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.shapiro.html

1

from scipy.stats import shapiro

经过Shapiro检验，发现p值较大，可以认为数据是来自于一个正态分布总体的

3.4.2 两独立样本均值检验¶

检验两个独立的样本对应总体的均值是否存在显著的差异

1

from scipy.stats import norm

1

n1 = norm(loc=0.3, scale=1.0)

1

x = np.linspace(-3,3,100)

1

from scipy.stats import ttest_ind #导入y检验模块

t = 2.1579961596973902p-value = 0.0321320234997519--------------------两总体均值存在显著差异

四、聚类分析—cluster¶

1

import scipy

随机产生20个4维点，disMat为各个点之间的距离矩阵

1

points=scipy.randn(20,4) #随机生成20个4维点

基于距离矩阵disMat进行层次聚类

1

# 进行层次聚类

从上图可以看到聚类的结果，例如我们考虑阈值为3，可以做一条y=3的水平线来辅助判断，可以看到图中9单独成一个类，17--11之间的点可形成一个类，剩余点也形成了一个类

1

cluster= sch.fcluster(Z, t=1) #t=1代表是阈值

层次聚类结果: [5 2 5 2 6 4 6 6 6 7 2 5 1 6 1 5 6 3 6 4]

五、插值处理——interpolate¶

5.1 一维插值¶

一般地通过多项式插值来进行插值，调用scipy.interpolate中的interp1d模块可以实现。通过kind参数可以调节插值的方法：

• nearest——最近邻插值
• zero——0阶插值
• linear——线性插值
• quadratic——二次插值
• cubic——三次插值
• 4,5。。。。——高阶多项式插值

1

x = np.arange(1,6)

[1 2 3 4 5][ 7.6883 19.5781 -9.014 9.7466 -7.3082]

1

from scipy.interpolate import interp1d

5.1.1 最近邻插值¶

1

plt.plot(X, ip_nearest(X), 'b*',label='最近邻插值点')

5.1.2 0阶插值¶

1

plt.plot(X, ip_zero(X), 'b*',label='0阶插值点')

5.1.3 线性插值¶

1

plt.plot(X, ip_linear(X), 'b*',label='线性插值点')

5.1.4 二次多项式插值¶

1

plt.plot(X, ip2(X), 'b*',label='二次多项式插值点')

5.1.5 三次多项式插值¶

1

plt.plot(X, ip3(X), 'b*',label='三次多项式插值点')

5.1.6 四次多项式插值¶

1

plt.plot(X, ip4(X), 'b*',label='四次多项式插值点')

5.2 二维差值¶

例如，相对下面的10$*$10的矩阵进行插值处理，扩张成一个100$*$100的矩阵，使得整个曲面更加光滑一些。

1

x1 = np.arange(1,11,1)

	1	2	3	4	5	6	7	8	9	10
1	-3.703201	-0.400318	1.580455	0.798721	6.071079	8.672773	6.820533	-0.707631	6.012882	0.380787
2	-1.145648	-3.759621	8.016974	4.090643	4.443695	0.446498	0.428496	3.707856	-17.316609	14.633428
3	-3.908394	-0.814071	1.315033	1.449469	-2.333758	4.574864	3.629699	-10.753554	-5.935987	-12.333878
4	-0.255300	-6.557950	-1.594488	-8.838964	5.145881	-0.901672	-5.951584	6.990636	11.796797	-7.305400
5	-9.913660	-12.400839	-6.915492	9.365875	-7.743060	-2.503458	3.416506	-6.089799	10.711856	18.799075
6	3.159555	-8.347973	-12.767721	0.614833	1.529912	-13.542717	1.809428	10.177388	7.766459	-6.101687
7	-8.817611	-0.471369	-3.194466	4.864112	0.914487	7.942217	-0.296812	-1.506718	-1.920830	4.727719
8	-15.056821	-1.742246	-5.434945	18.394791	0.529686	-10.661008	-7.973814	-0.112875	4.902229	-2.354906
9	-14.678328	-6.707785	5.278384	-7.699973	20.646028	-9.825168	-12.836315	-6.445505	-1.626839	-1.211078
10	36.102272	8.547623	-3.156888	-3.282717	8.829625	-15.125651	-15.816051	-2.890754	-5.458071	-17.222252

1

from scipy.interpolate import interp2d

1

interp = interp2d(x1, x2, y, kind='linear') #二维线性插值

1

X1 = np.linspace(x1.min(), x1.max(), 100)

1

# 可以看到插值后的 矩阵的维度为 100 * 100

1

from mpl_toolkits.mplot3d import Axes3D

5.3 径向基(RBF)插值¶

有些插值方法基于径向基函数，径向基函数的含义是点x处的函数值只依赖于x与某点c的距离。
常用的径向基函数(RBF)有Gaussian函数，Multiquadric函数和Inverse Multiquadric函数。

1

x = np.linspace(-3,3,100)

5.3.1 径向基函数介绍¶

1

plt.plot(x, np.exp(-1*x**2),linewidth=2,label=r'$y = e^{-x^2}$')

1

plt.plot(x, np.sqrt(1 + x **2),linewidth=2,label=r'$y = (1 + x^2)^{0.5}$')

1

plt.plot(x, 1. / np.sqrt(1 + x **2),linewidth=2,label=r'$y = \frac{1}{(1 + x^2)^{0.5}}$')

5.3.2 Rbf插值——以Guassian核为例¶

1

from scipy.interpolate.rbf import Rbf

1

x = np.arange(1,6)

[1 2 3 4 5][ -2.8767 4.0619 -15.9512 -16.7156 -7.4748]

1

X = np.arange(1,5,0.1)

5.3.3 高维RBF插值¶

1

x, y = np.mgrid[-np.pi/2:np.pi/2:5j, -np.pi/2:np.pi/2:5j]

x: [[-1.5708 -1.5708 -1.5708 -1.5708 -1.5708] [-0.7854 -0.7854 -0.7854 -0.7854 -0.7854] [ 0. 0. 0. 0. 0. ] [ 0.7854 0.7854 0.7854 0.7854 0.7854] [ 1.5708 1.5708 1.5708 1.5708 1.5708]]y: [[-1.5708 -0.7854 0. 0.7854 1.5708] [-1.5708 -0.7854 0. 0.7854 1.5708] [-1.5708 -0.7854 0. 0.7854 1.5708] [-1.5708 -0.7854 0. 0.7854 1.5708] [-1.5708 -0.7854 0. 0.7854 1.5708]]z: [[-0.6057 -0.1843 0. -0.1843 -0.6057] [-0.1843 0.444 0.7071 0.444 -0.1843] [ 0. 0.7071 1. 0.7071 0. ] [-0.1843 0.444 0.7071 0.444 -0.1843] [-0.6057 -0.1843 0. -0.1843 -0.6057]]

1

fig = plt.figure(figsize=(8,4))

原数据为三维空间上的一些散点，通过RBF插值可以形成一个光滑的曲面

1

zz = Rbf(x,y,z) #RBF插值

1

xx, yy = np.mgrid[-np.pi/2:np.pi/2:50j, -np.pi/2:np.pi/2:50j]

六、信号处理¶

6.1 序列线性趋势提取¶

例如对下面的x_series序列，希望提取其中的趋势项，并查看趋势项的形态

1

from scipy import signal

1

x1 = np.linspace(1,100,100)

对序列x_series做剔除趋势项处理，并将趋势项提取计算出来，查看趋势项的形态

1

after_trend = signal.detrend(x_series,type='linear') # 对序列x_series做剔除趋势项处理 type-linear指定趋势为线性的，也可设置为constant

6.2 傅里叶分析¶

傅里叶分析实际上是在频域考虑时域的问题，例如这里的将剔除某曲线(序列)中一些特定的频率成分，提取剩余部分成为滤波，再将频域中的滤波转化为时域当中显示的过程。
关于傅里叶分析可以参考：http://blog.jobbole.com/70549/
那么对于投资领域中股票的K线图，实际上是时域下的一条曲线，通过傅里叶分析将其转化到频域，再做频域下的处理，之后返回到时域上，可能会有出乎意料的效果。因为很多在时域上不好处理的操作，在频域上很容易实现。

1

from scipy import fftpack

利用傅里叶变换得到信号频谱，注意对象是剔除趋势项后的序列after_trend不是原序列

1

spec = np.abs(fftpack.fftshift(fftpack.rfft(after_trend)))

1

time = list(pd.date_range(start='2017-01-01', end='2017-04-10'))