python---scipy模块

一  简单介绍

SciPy是基于NumPy开发的高级模块，它提供了许多数学算法和函数的实现，用于解决科学计算中的一些标准问题。例如数值积分和微分方程求解，扩展的矩阵计算，最优化，概率分布和统计函数，甚至包括信号处理等。
作为标准科学计算程序库，SciPy类似于Matlab的工具箱，它是Python科学计算程序的核心包，它用于有效地计算NumPy矩阵，与NumPy矩阵协同工作。
SciPy库由一些特定功能的子模块构成，如下表所示：

模块	功能
cluster	矢量量化 / K-均值
constants	物理和数学常数
fftpack	傅里叶变换
integrate	积分程序
interpolate	插值
io	数据输入输出
linalg	线性代数程序
ndimage	n维图像包
odr	正交距离回归
optimize	优化
signal	信号处理
sparse	稀疏矩阵
spatial	空间数据结构和算法
special	任何特殊数学函数
stats	统计

以上子模块全依赖于NumPy且相互独立，导入NumPy和这些SciPy模块的标准方式如下，示例代码：

import numpy as np
from scipy import stats 

以上代码表示从SciPy模块中导入stats子模块，SciPy的其他子模块导入方式与之相同，限于机器学习研究领域及篇幅限制，本章将重点介绍linalg、optimize、interpolate及stats模块。

二 常用库的介绍

2.1 线性代数linalg模块

linalg是Linear Algebra的缩写，NumPy和SciPy都提供了线性代数函数库linalg，SciPy的线性代数库比NumPy更加全面。

（1）基本运算

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) 
0.0
>>> linalg.det(np.ones((3, 4)))        #无论行列式还是逆矩阵只适用于n阶矩阵的求解
Traceback (most recent call last):
...
ValueError: expected square matrix

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()函数用于比较两方阵所有对应元素值，如果完全相同返回真(True)，否则返回假(False)
True

以下计算奇异阵（行列式为0）的逆，其结果将会报错（LinAlgError），示例代码：

>>>arr = np.array([[3, 2], [6, 4]])
>>>linalg.inv(arr)
Traceback (most recent call last):
...
...LinAlgError: singular matrix

scipy.linalg.norm()函数计算方阵的范数，示例代码：

>>>A = np.matrix(np.random.random((2, 2)))
>>>A
>>>linalg.norm(A) #默认2范数
>>>linalg.norm(A, 1) #1范数
>>>linalg.norm(A, np.inf) #无穷范数

（2）解线性方程组

scipy.linalg.solve(A,b)和numpy.linalg.solve(A,b)可以用来解线性方程组Ax=b，即计算x=A^-1b。这里，A是mm的方形矩阵，x和b是长为m的向量。有时候A是固定的，需要对多组b进行求解，因此第二个参数也可以是mn的矩阵B。这样计算出来的X也是m*n的矩阵，相当于计算A^-1B。
在一些矩阵公式中经常会出现类似于A^-1B的运算，它们都可以用solve(A, B)计算，这要比直接逆矩阵然后做矩阵乘法更快捷一些，下面的程序比较solve()和逆矩阵的运算速度，示例代码：

>>> import numpy as np
>>> from scipy import linalg

>>> m, n = 500, 50
>>> A = np.random.rand(m, m)
>>> B = np.random.rand(m, n)
>>> X1 = linalg.solve(A, B)
>>> X2 = np.dot(linalg.inv(A), B)
>>> print(np.allclose(X1, X2))

>>> %timeit linalg.solve(A, B)
13.3 ms ± 834 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

>>> %timeit np.dot(linalg.inv(A), B)
22.4 ms ± 1.48 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

(3) 特征值和特征向量

n*n的矩阵A可以看作n维空间中的线性变换。若x为n维空间中的一个向量，那么A与x的矩阵乘积就是对x进行线性变换之后的向量。如果x是线性变换的特征向量，那么经过这个线性变换之后，得到的新向量仍然与原来的x保持在同一方向上，但其长度也许会改变。特征向量的长度在该线性变换下缩放的比例称为特征值。即特征向量x满足如下等式，λ的值可以是一个任意复数：Ax=λx。
下面以二维平面上的线性变换矩阵为例，演示特征值和特征向量的几何含义。通过linalg.eig(A)计算矩阵A的两个特征值evalues和特征向量evectors，在evectors中，每一列是一个特征向量。示例代码：

>>> A = np.array([[1, -0.3], [-0.1, 0.9]])
>>> evalues, evectors = linalg.eig(A)

2.2 拟合与求解optimize模块

SciPy的optimize模块提供了许多数值优化的算法，一些经典的优化算法包括线性回归、函数极值和根的求解以及确定两函数交点的坐标等。下面首先介绍简单的线性回归模型，然后逐渐深入解决非线性数据拟合问题。
（1）拟合 curve_fit()函数

线性回归有许多拟合数据的方法，我们将使用curve_fit()函数，它利用的是最小二乘算法。最小二乘算法是一种数学优化技术，在机器学习领域最有名和有效的算法之一。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。
以下示例中，我们首先从已知函数中生成一些带有噪声的数据，然后使用curve_fit()函数拟合这些噪声数据。示例中的已知函数我们使用一个简单的线性方程式，即f(x)=ax+b。示例代码：

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

#创建函数模型用来生成数据
def func(x, a, b):
       return a*x + b

#生成干净数据
x = np.linspace(0, 10, 100)
y = func(x, 1, 2)

#对原始数据添加噪声
yn = y + 0.9 * np.random.normal(size=len(x))

#使用curve_fit函数拟合噪声数据
popt, pcov = curve_fit(func, x, yn)

#输出给定函数模型func的最优参数
print(popt)

结果为：

[ 0.99734363  1.96064258]

如果有一个很好的拟合效果，popt返回的是给定模型的最优参数。我们可以使用pcov的值检测拟合的质量，其对角线元素值代表着每个参数的方差。

>>>print(pcov)
[[ 0.00105056 -0.00525282]
 [-0.00525282  0.03519569]]

通过以下代码绘制出了拟合曲线与实际曲线的差异，示例代码：

yfit = func(x,popt[0],popt[1]) 

plt.plot(x, y, color="green",label = "actual data")
plt.plot(x, yn, "o", label = "actual data with noise")
plt.plot(x, yfit,color="yellow", label = "fitting data")
plt.legend(loc = "best")
plt.show()

下面做进一步研究，我们可以通过最小二乘拟合高斯分布（Gaussian profile），一种非线性函数：α*exp(-(x-μ)²/2σ²)
在这里，α表示一个标量，μ是期望值，而σ是标准差。示例代码：

import numpy as np 
from scipy.optimize import curve_fit 
import matplotlib.pyplot as plt
#创建一个函数模型用来生成数据
def func(x, a, b, c):
        return (a*np.exp(-(x-b)**2/2*c**2))

#生成原始数据
x = np.linspace(0, 10, 100)
y = func(x, 1, 5, 2)

#对原始数据增加噪声
yn = y + 0.2*np.random.normal(size=len(x))

#使用curve_fit函数拟合噪声数据
popt, pcov = curve_fit(func, x, yn)

#popt返回最拟合给定的函数模型func的参数值，如popt[0]=a,popt[1]=b,popt[2]=3
print(popt)

结果为：

[-0.49627942  2.78765808 28.76127826]

通过以下代码绘制出了拟合曲线与实际曲线的差异，示例代码：

p0=[1.2,4,3] #初步猜测参数，如果没有，默认全为1，即[1,1,1]
popt, pcov = curve_fit(func, x, yn,p0=p0)

#popt返回最拟合给定的函数模型func的参数值，如popt[0]=a,popt[1]=b,popt[2]=3
print(popt)

yfit = func(x,popt[0],popt[1],popt[2])

plt.plot(x, y, color="green",label = "actual data")
plt.plot(x, yn, "o", label = "actual data with noise")
plt.plot(x, yfit, color="yellow", label = "fitting data")
plt.legend(loc = "best")
plt.show()

结果如下图所示：

 

 通过以上绘图，我们可以看出对高斯分布函数拟合的效果是可以接受的。
随着研究的深入，我们可以拟合一个多重高斯分布的一维数据集。现在将这个函数扩展为包含两个不同输入值的高斯分布函数。这是一个拟合线性光谱的经典实例，示例代码如下：

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
def func(x, a0, b0, c0, a1, b1, c1):
      return a0*np.exp(-(x - b0) ** 2/(2 * c0 ** 2)) + a1 * np.exp(-(x-b1) ** 2/(2 * c1 ** 2))

#生成原始数据
x = np.linspace(0, 20, 200)
y = func(x, 1, 3, 1, -2, 15, 0.5)

#对原始数据增加噪声
yn = y + 0.9 * np.random.normal(size=len(x))

#如果要拟合一个更加复杂的函数，提供一些估值假设对拟合效果更好
guesses = [1, 3, 1, 1, 15, 1]

#使用curve_fit函数拟合噪声数据
popt, pcov = curve_fit(func, x, yn, p0=guesses)

#popt返回最拟合给定的函数模型func的参数值，如popt[0]=a,popt[1]=b,popt[2]=3
print(popt)

yfit = func(x,popt[0],popt[1],popt[2],popt[3],popt[4],popt[5])

plt.plot(x, y, color="green",label = "actual data")
plt.plot(x, yn, "o", label = "actual data with noise")
plt.plot(x, yfit, color="yellow", label = "fitting data")
plt.legend(loc = "best")
plt.show()

结果如下图所示：

（2）最小二乘拟合leastsq()函数

假设有一组实验数据(x[i], y[i])，我们知道它们之间的函数关系:y = f(x)，通过这些已知信息，需要确定函数中的一些参数项。例如，如果f是一个线型函数f(x) = k*x+b，那么参数k和b就是我们需要确定的值。如果将这些参数用 p 表示的话，那么我们就是要找到一组 p 值使得如下公式中的S函数最小：

这种算法被称之为最小二乘拟合(Least-square fitting)。optimize模块提供了实现最小二乘拟合算法的函数leastsq()，leastsq是least square的简写，即最小二乘法。下面是用leastsq()对线性函数进行拟合的程序，示例代码：

import matplotlib.pylab as plt
import numpy as np
from scipy import optimize    # 从scipy库引入optimize模块

X = np.array([ 8.19, 2.72, 6.39, 8.71, 4.7, 2.66, 3.78 ])
Y = np.array([ 7.01, 2.78, 6.47, 6.71, 4.1, 4.23, 4.05 ])

def residuals(p):
        #计算以p为参数的直线和原始数据之间的误差
        k, b = p
        return Y-(k*X+b)

# leastsq()使得residuals()的输出数组的平方和最小，参数的初始值为[1, 0]
r = optimize.leastsq(residuals, [1,0])
k, b = r[0]
print("k=", k, "b=", b)

结果为：

k = 0.613495349193  b = 1.79409254326

可以通过通过绘图对比真实数据和拟合数据的误差，示例代码；

plt.plot(X, Y, "o", label = "actual data")
plt.plot(X, k*X+b, label = "fitting data")
plt.legend(loc = "best")
plt.show()

结果为：

绘图中的圆点表示真实数据点，实线表示拟合曲线，由此看出拟合参数得到的函数和真实数据大体一致。接下来，用leastsq()对正弦波数据进行拟合，示例代码：

import numpy as np
from scipy.optimize import leastsq   # 从scipy库的optimize模块引入leastsq函数
import matplotlib.pyplot as plt    # 引入绘图模块pylab，并重命名为pl

def func(x, p):
    """
    数据拟合所用的函数: A*sin(2*pi*k*x + theta)
    """
    A, k, theta = p
    return A*np.sin(2*np.pi*k*x+theta)   

def residuals(p, y, x):
    """
    实验数据x, y和拟合函数之间的差，p为拟合需要找到的系数
    """
    return y - func(x, p) 

x = np.linspace(0, -2*np.pi, 100)
A, k, theta = 10, 0.34, np.pi/6   # 真实数据的函数参数
y0 = func(x, [A, k, theta])   # 真实数据

y1 = y0 + 2 * np.random.randn(len(x))   # 加入噪声之后的实验数据，噪声是服从标准正态分布的随机量    

p0 = [7, 0.2, 0]   # 第一次猜测的函数拟合参数

# 调用leastsq进行数据拟合
# residuals为计算误差的函数
# p0为拟合参数的初始值
# args为需要拟合的实验数据
plsq = leastsq(residuals, p0, args=(y1, x))

print ("actual parameter:", [A, k, theta]) # 真实参数
print ("fitting parameter", plsq[0]) # 实验数据拟合后的参数

plt.plot(x, y0, label="actual data") # 绘制真实数据
plt.plot(x, y1, label="experimental data with noise")  # 带噪声的实验数据
plt.plot(x, func(x, plsq[0]), label="fitting data")    # 拟合数据
plt.legend()
plt.show()

这个例子中我们要拟合的函数是一个正弦波函数，它有三个参数 A, k, theta ，分别对应振幅、频率、相角。假设我们的实验数据是一组包含噪声的数据 x, y1，其中y1是在真实数据y0的基础上加入噪声的到了。通过leastsq函数对带噪声的实验数据x, y1进行数据拟合，可以找到x和真实数据y0之间的正弦关系的三个参数： A, k, theta。下面是程序的输出：

>>>actual parameter: [10, 0.34, 0.5235987755982988]
>>>fitting parameter [ 10.12646889   0.33767587   0.48944317]

我们看到拟合参数虽然和真实参数完全不同，但是由于正弦函数具有周期性，实际上拟合参数得到的函数和真实参数对应的函数是一致的。
（3）标量函数极值求解fmin()函数
首先定义以下函数，然后绘制它，示例代码：

import numpy as np
from scipy import optimize
import matplotlib.pyplot as plt
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算法是做这个的好方法，BFGS算法被认为是数值效果最好的拟牛顿法，是由Broyden，Fletcher，Goldfarb，Shanno四个人分别提出的，故称为BFGS校正。具体算法思想及解释请查阅相关资料，这里直接通过optimize.fmin_bfgs()函数求解最小值，示例代码：

>>> optimize.fmin_bfgs(f, 0)
Optimization terminated successfully.
         Current function value: -7.945823
         Iterations: 5
         Function evaluations: 24
         Gradient evaluations: 8
array([-1.30644003])

这个方法一个可能的问题在于，如果函数有局部最小值，算法会因初始点不同找到这些局部最小而不是全局最小，示例代码：

>>> optimize.fmin_bfgs(f, 3, disp=0)#disp是布尔型数据，如果为1，打印收敛消息
array([ 3.83746663])

如果我们不知道全局最小值的邻近值来选定初始点，我们需要借助于耗费资源些的全局优化。为了找到全局最小点，最简单的算法是蛮力算法，该算法求出给定格点的每个函数值。示例代码：

>>>grid = (-10, 10, 0.1)
>>>xmin_global = optimize.brute(f, (grid, ))
>>>xmin_global
array([-1.30641113])

对于大点的格点，scipy.optimize.brute()变得非常慢。scipy.optimize.anneal()提供了使用模拟退火的替代函数。对已知的不同类别全局优化问题存在更有效率的算法，但这已经超出scipy的范围。为了找到局部最小，我们把变量限制在(0,10)之间，使用scipy.optimize.fminbound()，示例代码：

>>> xmin_local = optimize.fminbound(f, 0, 10)
>>> xmin_local
3.8374671...

下面的程序通过求解卷积的逆运算演示fmin的功能。对于一个离散线性时不变系统h, 如果输入是x，那么其输出y可以用x和h的卷积表示：

现在的问题是如果已知系统的输入x和输出y，如何计算系统的传递函数h；或者如果已知系统的传递函数h和系统的输出y，如何计算系统的输入x。这种运算被称为反卷积运算，是十分困难的，特别是在实际的运用中，测量系统的输出总是存在误差的。下面用fmin计算反卷积，这种方法只能用在很小规模的数列之上，因此没有很大的实用价值，不过用来评价fmin函数的性能还是不错的。示例代码：

import scipy.optimize as opt 
import numpy as np 

def test_fmin_convolve(fminfunc, x, h, y, yn, x0): 
    """
    x (*) h = y, (*)表示卷积
    yn为在y的基础上添加一些干扰噪声的结果
    x0为求解x的初始值
    """
    def convolve_func(h):
        """
        计算 yn - x (*) h 的power
        fmin将通过计算使得此power最小
        """ 
        return np.sum((yn - np.convolve(x, h))**2) 

    # 调用fmin函数，以x0为初始值
    h0 = fminfunc(convolve_func, x0) 

    print fminfunc.__name__ 
    print "---------------------" 
    # 输出 x (*) h0 和 y 之间的相对误差
    print "error of y:", np.sum((np.convolve(x, h0)-y)**2)/np.sum(y**2) 
    # 输出 h0 和 h 之间的相对误差
    print "error of h:", np.sum((h0-h)**2)/np.sum(h**2) 
    print 

def test_n(m, n, nscale): 
    """
    随机产生x, h, y, yn, x0等数列，调用各种fmin函数求解b
    m为x的长度, n为h的长度, nscale为干扰的强度
    """
    x = np.random.rand(m) 
    h = np.random.rand(n) 
    y = np.convolve(x, h) 
    yn = y + np.random.rand(len(y)) * nscale
    x0 = np.random.rand(n) 

    test_fmin_convolve(opt.fmin, x, h, y, yn, x0) 
    test_fmin_convolve(opt.fmin_powell, x, h, y, yn, x0) 
    test_fmin_convolve(opt.fmin_cg, x, h, y, yn, x0)
    test_fmin_convolve(opt.fmin_bfgs, x, h, y, yn, x0)

test_n(200, 20, 0.1)

运行结果为：

fmin
---------------------
error of y: 0.000360456186137
error of h: 0.0122264525455
Optimization terminated successfully.
         Current function value: 0.207509
         Iterations: 96
         Function evaluations: 17400
fmin_powell
---------------------
error of y: 0.000129249083036
error of h: 0.000300953639205
Optimization terminated successfully.
         Current function value: 0.207291
         Iterations: 20
         Function evaluations: 880
         Gradient evaluations: 40
fmin_cg
---------------------
error of y: 0.000129697740414
error of h: 0.000292820536053
Optimization terminated successfully.
         Current function value: 0.207291
         Iterations: 31
         Function evaluations: 946
         Gradient evaluations: 43
fmin_bfgs
---------------------
error of y: 0.000129697643272
error of h: 0.000292817401206

（4）函数求解fsolve()

optimize库中的fsolve函数可以用来对非线性方程组进行求解，其基本调用形式是

                                                                              fsolve(func, x0)

• func是用于定义需求解的非线性方程组的函数文件名
• x0为未知数矢量的初始值。

首先通过一个简单的示例，利用fsolve()函数求解当线性函数为0时，x的值，示例代码：

import matplotlib.pyplot as plt
from scipy.optimize import fsolve
import numpy as np

line = lambda x:x+3

solution = fsolve(line, -2)
print(solution)

结果为：

[-3，]

通过以下绘图函数可以看出当函数等于0时，x轴的坐标值为-3，示例代码：

x = np.linspace(-5.0, 0, 100)
plt.plot(x,line(x), color="green",label = "function")
plt.plot(solution,line(solution), "o", label = "root")
plt.legend(loc = "best")
plt.show()

结果为：

下面我们通过一个简单的示例介绍一下两个方程交点的求解方法，示例代码：

from scipy.optimize import fsolve
import numpy as np
import matplotlib.pyplot as plt
# 定义解函数
def findIntersection(func1, func2, x0):
        return fsolve(lambda x: func1(x)-func2(x),x0)

# 定义两方程
funky = lambda x : np.cos(x / 5) * np.sin(x / 2)
line = lambda x : 0.01 * x - 0.5

# 定义两方程交点的取值范围
x = np.linspace(0, 45, 1000)
result = findIntersection(funky, line, [15, 20, 30, 35, 40, 45])

# 输出结果
print(result, line(result))


plt.plot(x,funky(x), color="green",label = "funky func")
plt.plot(x,line(x), color="yellow",label = "line func")
plt.plot(result,line(result), "o", label = "intersection")
plt.legend(loc = "best")
plt.show()

结果为：

如果要对如下方程组进行求解的话：

• f1(u1,u2,u3) = 0
• f2(u1,u2,u3) = 0
• f3(u1,u2,u3) = 0

那么func可以如下定义：

def func(x):
    u1,u2,u3 = x
    return [f1(u1,u2,u3), f2(u1,u2,u3), f3(u1,u2,u3)]

下面是一个实际的例子，求解如下方程组的解：

• 5*x1 + 3 = 0
• 4*x0*x0 - 2*sin(x1*x2) = 0
• x1*x2 - 1.5 = 0

示例代码：

from scipy.optimize import fsolve
from math import sin,cos

def f(x):
    x0 = float(x[0])
    x1 = float(x[1])
    x2 = float(x[2])
    return [
        5*x1+3,
        4*x0*x0 - 2*sin(x1*x2),
        x1*x2 - 1.5
    ]
result = fsolve(f, [1,1,1])
print (result)

结果为：

[-0.70622057 -0.6        -2.5       ]

2.3 插值interpolate模块

插值是离散函数逼近的重要方法，利用它可通过函数在有限个点处的取值状况，估算出函数在其他点处的近似值。与拟合不同的是，要求曲线通过所有的已知数据。SciPy的interpolate模块提供了许多对数据进行插值运算的函数，范围涵盖简单的一维插值到复杂多维插值求解。当样本数据变化归因于一个独立的变量时，就使用一维插值；反之样本数据归因于多个独立变量时，使用多维插值。
计算插值有两种基本的方法，1、对一个完整的数据集去拟合一个函数；2、对数据集的不同部分拟合出不同的函数，而函数之间的曲线平滑对接。第二种方法又叫做仿样内插法，当数据拟合函数形式非常复杂时，这是一种非常强大的工具。我们首先介绍怎样对简单函数进行一维插值运算，然后进一步深入比较复杂的多维插值运算。

（1）一维插值

一维数据的插值运算可以通过函数interp1d()完成。其调用形式如下，它实际上不是函数而是一个类：

interp1d(x, y, kind='linear', ...)

其中，x和y参数是一系列已知的数据点，kind参数是插值类型，可以是字符串或整数，它给出插值的B样条曲线的阶数，候选值及作用下表所示：

候选值	作用
‘zero’ 、'nearest'	阶梯插值，相当于0阶B样条曲线
‘slinear’ 、'linear'	线性插值，用一条直线连接所有的取样点，相当于一阶B样条曲线
‘quadratic’ 、'cubic'	二阶和三阶B样条曲线，更高阶的曲线可以直接使用整数值指定

下面的程序演示了通过不同的 kind参数（linear和quadratic），对一个正弦函数进行插值运算。示例代码：

import numpy as np
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt

#创建待插值的数据
x = np.linspace(0, 10*np.pi, 20)
y = np.cos(x)

# 分别用linear和quadratic插值
fl = interp1d(x, y, kind='linear')
fq = interp1d(x, y, kind='quadratic')

#设置x的最大值和最小值以防止插值数据越界
xint = np.linspace(x.min(), x.max(), 1000)
yintl = fl(xint)
yintq = fq(xint)


plt.plot(xint,fl(xint), color="green", label = "Linear")
plt.plot(xint,fq(xint), color="yellow", label ="Quadratic")
plt.legend(loc = "best")
plt.show()

结果如下图所示：

（2）噪声数据插值

 我们可以通过interpolate模块中UnivariateSpline()函数对含有噪声的数据进行插值运算，示例代码：

import numpy as np
from scipy.interpolate import UnivariateSpline
import matplotlib.pyplot as plt

# 通过人工方式添加噪声数据
sample = 30
x = np.linspace(1, 10*np.pi, sample)
y = np.cos(x) + np.log10(x) + np.random.randn(sample)/10

# 插值，参数s为smoothing factor 
f = UnivariateSpline(x, y, s=1)
xint = np.linspace(x.min(), x.max(), 1000)
yint = f(xint)

plt.plot(xint,f(xint), color="green", label = "Interpolation")
plt.plot(x, y, color="yellow", label ="Original")
plt.legend(loc = "best")
plt.show()

需要说明的是：在UnivariateSpline()函数中，参数s是平滑向量参数，被用来拟合还有噪声的数据。如果参数s=0，将忽略噪声对所有点进行插值运算。结果如下图所示：

 （3）多维插值

 多维插值主要用于重构图片，interpolate模块中的griddata()函数有很强大的处理多维散列取样点进行插值运算的能力。其调用形式如下：

griddata(points, values, xi, method='linear', fill_value=nan)

其中points表示K维空间中的坐标，它可以是形状为(N,k)的数组，也可以是一个有k个数组的序列，N为数据的点数。values是points中每个点对应的值。xi是需要进行插值运算的坐标，其形状为(M,k)。method参数有三个选项：'nearest'、 ‘linear’、 'cubic'，分别对应0阶、1阶以及3阶插值。以下示例利用1000个随机散列点对1000x1000像素的图片进行重构，示例代码：

import numpy as np
from scipy.interpolate import griddata#定义一个函数
def ripple(x,y):
    return np.sqrt(x**2 + y**2) + np.sin(x**2 + y**2)

# 生成grid数据，复数定义了生成grid数据的step，若无该复数则step为5 
grid_x, grid_y = np.mgrid[0:5:1000j, 0:5:1000j] 

# 生成待插值的样本数据 
points = np.random.rand(1000,2) 

value = ripple(points[:,0]*5,points[:,1]*5) 

# 用nearest方法插值
grid_z0 = griddata(points*5,value, (grid_x,grid_y),method='nearest')

我们还可以使用interpolate模块的SmoothBivariateSpline类进行多元仿样插值运算，对图片进行重构。示例代码：

import numpy as np
from scipy.interpolate import SmoothBivariateSpline as SBS

#定义一个函数
def ripple(x,y):
    return np.sqrt(x**2 + y**2) + np.sin(x**2 + y**2)

# 生成grid数据，复数定义了生成grid数据的step，若无该复数则step为5 
grid_x, grid_y = np.mgrid[0:5:1000j, 0:5:1000j] 

# 生成待插值的样本数据 
points = np.random.rand(1000,2)
 
value = ripple(points[:,0]*5,points[:,1]*5) 

# 用nearest方法插值
fit = SBS(points[:,0]*5, points[:,1]*5, value, s=0.01, kx=4, ky=4)
interp = fit(np.linspace(0, 5, 1000), np.linspace(0, 5, 1000))

我们得到了一个与上个示例同样的结果。整体上SmoothBivariateSpline函数的表现略好于griddata函数。
通过反复测试，尽管SmoothBivariateSpline表现略好，但其对给定的样本数据非常敏感，就可能导致忽略一些显著特征。而griddata函数有很强的鲁棒性，不管给定的数据样本，能够合理的进行插值运算。

2.4 统计stats模块
NumPy库已经提供了一些基本的统计函数，如求期望、方差、中位数、最大值和最小值等。示例代码：

import numpy as np

#构建一个1000个随机变量的数组
x = np.random.randn(1000)

#对数组元素的值进行统计
mean = x.mean()
std = x.std()
var = x.var()

print(mean,std,var)

结果为：

(0.02877273942510088, 0.97623362287515114, 0.95303208643194282)

mean是期望值，std是标准差，var是方差，使用numpy.array对象已有的方法获得统计指标快速有效，而SciPy库则提供了更高级的统计工具，它的Stats模块包含了多种概率分布的随机变量（随机变量是指概率论中的概念，不是Python中的变量），其中随机变量又分为连续和离散两种。所有的连续随机变量都是rv_continuous的派生类的对象，而所有的离散随机变量都是rv_discrete的派生类的对象。

（1）连续概率分布
SciPy的stats模块提供了大约80种连续随机变量和10多种离散分布变量，这些分布都依赖于numpy.random函数。可以通过如下语句获得stats模块中所有的连续随机变量，示例代码：

from scipy import stats
[k for k, v in stats.__dict__.items() if isinstance(v, stats.rv_continuous)]

运行结果为：

['ksone', 'kstwobign', 'norm', 'alpha', 'anglit', 'arcsine', 'beta', 'betaprime', 'bradford', 'burr', 'burr12', 'fisk', 'cauchy', 'chi', 'chi2', 'cosine', 'dgamma', 'dweibull', 'expon', 'exponnorm', 'exponweib', 'exponpow', 'fatiguelife', 'foldcauchy', 'f', 'foldnorm', 'frechet_r', 'weibull_min', 'frechet_l', 'weibull_max', 'genlogistic', 'genpareto', 'genexpon', 'genextreme', 'gamma', 'erlang', 'gengamma', 'genhalflogistic', 'gompertz', 'gumbel_r', 'gumbel_l', 'halfcauchy', 'halflogistic', 'halfnorm', 'hypsecant', 'gausshyper', 'invgamma', 'invgauss', 'invweibull', 'johnsonsb', 'johnsonsu', 'laplace', 'levy', 'levy_l', 'levy_stable', 'logistic', 'loggamma', 'loglaplace', 'lognorm', 'gilbrat', 'maxwell', 'mielke', 'kappa4', 'kappa3', 'nakagami', 'ncx2', 'ncf', 't', 'nct', 'pareto', 'lomax', 'pearson3', 'powerlaw', 'powerlognorm', 'powernorm', 'rdist', 'rayleigh', 'reciprocal', 'rice', 'recipinvgauss', 'semicircular', 'skewnorm', 'trapz', 'triang', 'truncexpon', 'truncnorm', 'tukeylambda', 'uniform', 'vonmises', 'vonmises_line', 'wald', 'wrapcauchy', 'gennorm', 'halfgennorm']

连续随机变量对象主要使用如下方法，下表所示：

方法名	全称	功能
rvs	Random Variates of given type	对随机变量进行随机取值，通过size参数指定输出数组的大小
pdf	Probability Density Function	随机变量的概率密度函数
cdf	Cumulative Distribution Function	随机变量的累积分布函数，它是概率密度函数的积分
sf	Survival function	随机变量的生存函数，它的值是1-cdf(t)
ppf	Percent point function	累积分布函数的反函数
stats	statistics	计算随机变量的期望值和方差
fit	fit	对一组随机取样进行拟合，找出最适合取样数据的概率密度函数的系数

下面以标准正态分布（函数表示f(x)=(1/√2π)exp(-x^2/2)）为例，简单介绍随机变量的用法。示例代码：

from scipy import stats
# 设置正态分布参数，其中loc是期望值参数，scale是标准差参数
X = stats.norm(loc=1.0, scale=2.0)

# 计算随机变量的期望值和方差
print(X.stats())

结果为：

(array(1.0), array(4.0))

以上代码说明，norm可以像函数一样调用，通过loc和scale参数可以指定随机变量的偏移和缩放参数。对于正态分布的随机变量来说，这两个参数相当于指定其期望值和标准差，标准差是方差的算术平方根。X的stats()方法，可以计算随机变量X分布的特征值，如期望值和方差。
此外，通过调用随机变量X的rvs()方法，可以得到包含一万次随机取样值的数组x，然后调用NumPy的mean()和var()计算此数组的均值和方差，其结果符合随机变量X的特性，示例代码：

#对随机变量取10000个值
x = X.rvs(size=10000)
print(np.mean(x), np.var(x))

结果为：

(1.0287787687588861, 3.9944276709242805)

使用fit()方法对随机取样序列x进行拟合，它返回的是与随机取样值最吻合的随机变量参数，示例代码：

#输出随机序列的期望值和标准差
print(stats.norm.fit(x))

结果为：

(1.0287787687588861, 1.998606432223283)

在下面的例子中，计算取样值x的直方图统计以及累计分布，并与随机变量的概率密度函数和累积分布函数进行比较。示例代码：

pdf, t = np.histogram(x, bins=100, normed=True)
t = (t[:-1]+t[1:])*0.5
cdf = np.cumsum(pdf) * (t[1] - t[0])
p_error = pdf - X.pdf(t)
c_error = cdf - X.cdf(t)
print("max pdf error: {}, max cdf error: {}".format(np.abs(p_error).max(), np.abs(c_error).max()))

运行结果如下所示：

max pdf error: 0.0208405611169, max cdf error: 0.0126874590568

通过绘图的方式查看概率密度函数求得的理论值（theory value）和直方图统计值（statistic value），可以看出二者是一致的，示例代码：

import pylab as pl
pl.plot(t, pdf, color="green", label = "statistic value")
pl.plot(t, X.pdf(t), color="yellow", label ="theory value")
pl.legend(loc = "best")
pl.show()

结果见下图所示：

也可以用同样的方式显示随机变量X的累积分布和数组pdf的累加结果，示例代码：

import pylab as pl
pl.plot(t, cdf, color="green", label = "statistic value")
pl.plot(t, X.cdf(t), color="yellow", label ="theory value")
pl.legend(loc = "best")
pl.show()

结果为：

（2）离散概率分布

# 数组x保存骰子的所有可能值，数组p保存每个值出现的概率
x = range(1, 7)
p = (0.4, 0.2, 0.1, 0.1, 0.1, 0.1)

# 创建表示这个骰子的随机变量dice，调用其rvs()方法投掷此骰子20次，获得符合概率p的随机数
dice = stats.rv_discrete(values=(x, p))
print(dice.rvs(size=20))

运行结果：

array([3, 6, 4, 5, 5, 2, 1, 3, 3, 1, 1, 3, 1, 5, 1, 3, 4, 1, 2, 2])

除了自定义离散概率分布，我们也可以利用stats模块里的函数定义各种分布。下面以生成几何分布为例，其函数是geom()，示例代码：

import numpy as np
from scipy.stats import geom

# 设置几何分布的参数
p = 0.5
dist = geom(p)

# 设置样本区间  
x = np.linspace(0, 5, 1000)  

# 得到几何分布的 PMF 和CDF  
pmf = dist.pmf(x) 
cdf = dist.cdf(x)  

# 生成500个随机数  
sample = dist.rvs(500)

（3）描述与检验函数

 SciPy中有超过60种统计函数。stats模块包括了诸如kstest 和normaltest等样本测试函数，用来检测样本是否服从某种分布。在使用这些工具前，要对数据有较好的理解，否则可能会误读它们的结果。样本分布检验为例，示例代码：

import numpy as np 
from scipy import stats 

# 生成包括100个服从正态分布的随机数样本
sample = np.random.randn(100) 

# 用normaltest检验原假设
out = stats.normaltest(sample) 
print('normaltest output') 
print('Z-score = ' + str(out[0])) 
print('P-value = ' + str(out[1])) 

# kstest 是检验拟合度的Kolmogorov-Smirnov检验，这里针对正态分布进行检验
# D是KS统计量的值，越接近0越好
out = stats.kstest(sample, 'norm') 
print('
kstest output for the Normal distribution') 
print('D = ' + str(out[0])) 
print('P-value = ' + str(out[1])) 

# 类似地可以针对其他分布进行检验，例如Wald分布
out = stats.kstest(sample, 'wald') 
print('
kstest output for the Wald distribution') 
print('D = ' + str(out[0])) 
print('P-value = ' + str(out[1]))

SciPy的stats模块中还提供了一些描述函数，如几何平均（gmean）、偏度（skew）、样本频数（itemfreq）等。示例代码

import numpy as np 
from scipy import stats 

# 生成包括100个服从正态分布的随机数样本
sample = np.random.randn(100) 

# 调和平均数，样本值须大于0 
out = stats.hmean(sample[sample > 0]) 
print('Harmonic mean = ' + str(out)) 

# 计算-1到1之间样本的均值
out = stats.tmean(sample, limits=(-1, 1)) 
print('
Trimmed mean = ' + str(out)) 

# 计算样本偏度
out = stats.skew(sample) 
print('
Skewness = ' + str(out)) 

# 函数describe可以一次给出样本的多种描述统计结果
out = stats.describe(sample) 
print('
Size = ' + str(out[0])) 
print('Min = ' + str(out[1][0])) 
print('Max = ' + str(out[1][1])) 
print('Mean = ' + str(out[2])) 
print('Variance = ' + str(out[3])) 
print('Skewness = ' + str(out[4])) 
print('Kurtosis = ' + str(out[5]))

参考：简书