Scipy在NumPy的基础上增加了众多的数学计算、科学计算以及工程计算中常用的模块。接下来会对其慢慢介绍
一、常数和函数
(一)常数
Scipy的constants
模块包含了众多的物理常数
from scipy import constants as C
print(C.c) # 真空中的光速
print(C.h) # 普朗克常数
在字典physical_constants
中,以物理常量名为键,对应值是一个含有三个元素的元组,分别对应常数值、单位和误差
# 查看电子的质量
from scipy import constants as C
print(C.physical_constants["electron mass"])
除了物理常数之外,constans
模块还包括许多单位信息,它们是1单位的量转换成标准单位时的数值
print(C.mile) # 1英里等于多少米
print('')
print(C.inch) # 1英寸等于多少米
print('')
print(C.gram) # 1克等于多少千克
print('')
print(C.pound) # 1磅等于多少千克
(二)函数
special
模块是一个非常完整的函数库,其中包含了基本数学函数、特殊数学函数以及NumPy中出现的所有函数
伽马函数
Γ
\Gamma
Γ是概率统计中常出现的一个特殊函数:
Γ ( z ) = ∫ 0 ∞ t z − 1 e − t d t \Gamma(z)=\int_{0}^{\infty}t^{z-1}e^{-t}dt Γ(z)=∫0∞tz−1e−tdt
在SciPy中可以通过special模块中的gamma()进行计算
import scipy.special as S
print(S.gamma(4))
print('')
print(S.gamma(0.5))
print('')
print(S.gamma(1+1j)) # gamma函数支持复数
print('')
print(S.gamma(1000))
special模块中的某些函数并不是数学意义上的特殊函数,例如logp(x)
计算log(1+x)
的值,这是由于浮点数的精度有限导致的,当使用loglp()
时,就可以解决这个问题
print(1 + 1e-20)
print('')
print(np.log(1+1e-20))
print('')
print(S.log1p(1e-20))
二、optimize
(一)非线性方程组求解
fsolve()
可以对非线性方程组进行求解:
fsolve(func, x0)
其中func是计算方程组误差的函数,它的参数x是一个数组,其值为方程组的一组可能的解。func返回将x代入方程组之后得到的每个方程的误差,x0为未知数的一组初始值。
from math import sin, cos
from scipy import optimize
def f(x):
x0, x1, x2 = x.tolist()
return [
5*x1+3,
4*x0*x0 - 2*sin(x1*x2),
x1*x2 - 1.5
]
# f计算方程组的误差,[1,1,1]是未知数的初始值
result = optimize.fsolve(f, [1,1,1])
print(result)
print(f(result))
该过程实现了对方程组:
[
5
x
+
3
=
0
4
x
0
2
−
2
s
i
n
(
x
1
x
2
)
=
0
x
1
x
2
−
1.5
=
0
]
\begin{bmatrix} 5x+3=0\\ 4x_0^2-2sin(x_1x_2)=0\\ x_1x_2-1.5=0 \end{bmatrix}
⎣⎡5x+3=04x02−2sin(x1x2)=0x1x2−1.5=0⎦⎤的求解
在对方程组进行求解时,fsolve()
会自动计算方程组在某点对各个未知数变量的偏导数,这些偏导数组成一个二维数组,也就是雅可比矩阵。
如果方程组中的未知数很多,而与每个方程有关联的未知数较少,即雅可比矩阵比较稀疏时,将计算雅可比矩阵的函数作为参数传给fsolve()
,能大幅度提高计算速度。
(二)最小二乘拟合
假如有一组数据(xi, yi)且假设它们满足某种函数关系yi = f(xi)
。
若f()=kx+b
,为了确定k和b的值,我们可以用最小二乘拟合来确定
如果用p表示函数中需要确定的参数,则目标是找到一组p使得函数s的值最小:
S
(
P
)
=
∑
i
=
1
m
[
y
i
−
f
(
x
i
,
p
)
]
2
S(P)=\sum_{i=1}^{m}[y_i-f(x_i,p)]^2
S(P)=∑i=1m[yi−f(xi,p)]2
这种算法就称作最小二乘拟合(least_square fitting)。
在optimize模块中,可以使用leastsq()
对数据进行最小二乘拟合计算
import numpy as np
from scipy import 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 ={}\n".format(k), "b ={}\n".format(b))
(三)计算函数局域最小值
optimize库还提供了许多求函数最小值的方法:Nelder-Mead、Powell、CG、BFGS、Newton-CG、L-BFGS-B等。
为了提高运算速度和精度,有些算法带有一个fprime参数,用来计算目标函数f()对各个自变量的偏导数的函数
在这取f(x, y)为:
f
(
x
,
y
)
=
(
1
−
x
)
2
+
100
(
y
−
x
2
)
2
f(x,y)=(1-x)^2+100(y-x^2)^2
f(x,y)=(1−x)2+100(y−x2)2
则f(x, y)对x和y的偏导函数为:
∂ f ∂ x = − 2 + 2 x − 400 x ( y − x 2 ) \frac{\partial f}{\partial x}=-2+2x-400x(y-x^2) ∂x∂f=−2+2x−400x(y−x2)
∂ f ∂ y = 200 y − 200 x 2 \frac{\partial f}{\partial y}=200y-200x^2 ∂y∂f=200y−200x2
而Newton-CG算法需要计算海森矩阵,它是一个由自变量为向量的实值函数的二阶偏导数构成的方块矩阵
则对于该f(x, y)来说,海森矩阵为一个二阶矩阵:
[ 2 ( 600 x 2 − 200 y + 1 ) − 400 x − 400 x 200 ] \begin{bmatrix} 2(600x^2-200y+1)&-400x \\ -400x& 200 \end{bmatrix} [2(600x2−200y+1)−400x−400x200]
最后用各种最小值优化算法计算f(x, y)的最小值,根据结果可以看出哪些算法需要较长的收敛时间,哪些算法能更快找到最小点
import numpy as np
from scipy import optimize
def target_function(x, y):
return (1-x)**2 + 100*(y-x**2)**2
class TargetFunction(object):
def __init__(self):
self.f_points = []
self.fprime_points = []
self.fhess_points = []
def f(self, p):
x, y = p.tolist()
z = target_function(x, y)
self.f_points.append((x, y))
return z
def fprime(self, p):
x, y = p.tolist()
self.fprime_points.append((x, y))
dx = -2 + 2*x - 400*x*(y - x**2)
dy = 200*y - 200*x**2
return np.array([dx, dy])
def fhess(self, p):
x, y = p.tolist()
self.fhess_points.append((x, y))
return np.array([[2*(600*x**2 - 200*y + 1), -400*x],
[-400*x, 200]])
def fmin_demo(method):
target = TargetFunction()
init_point =(-1, -1)
res = optimize.minimize(target.f, init_point,
method=method,
jac=target.fprime,
hess=target.fhess)
return res, [np.array(points) for points in
(target.f_points, target.fprime_points, target.fhess_points)]
methods = ("Nelder-Mead", "Powell", "CG", "BFGS", "Newton-CG", "L-BFGS-B")
for method in methods:
res, (f_points, fprime_points, fhess_points) = fmin_demo(method)
print("{:12s}: min={:12g}, f count={:3d}, fprime count={:3d}, fhess count={:3d}".format(
method, float(res["fun"]), len(f_points), len(fprime_points), len(fhess_points)))
(四)计算函数全域最小值
前面介绍的集中最小值优化算法只能计算局域的最小值,optimize库还提供了集中能进行全局优化的算法
下面一正弦波拟合为例,演示全局优化函数的算法
def func(x, p):
A, k, theta = p
return A*np.sin(2*np.pi*k*x+theta)
def func_error(p, y, x):
return np.sum((y - func(x, p))**2)
x = np.linspace(0, 2*np.pi, 100)
A, k, theta = 10, 0.34, np.pi/6
y0 = func(x, [A, k, theta])
np.random.seed(0)
y1 = y0 + 2 * np.random.randn(len(x))
result = optimize.basinhopping(func_error, (1, 1, 1),
niter = 10,
minimizer_kwargs={"method":"L-BFGS-B",
"args":(y1, x)})
print(result.x)
由于正弦函数的周期性,其拟合曲线是和原始数据重合的
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
plt.rcParams['font.sans-serif']=['SimHei'] # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus']=False # 用来正常显示负号
#%figonly=用`basinhopping()`拟合正弦波
plt.plot(x, y1, "o", label=u"带噪声的实验数据")
plt.plot(x, y0, label=u"真实数据")
plt.plot(x, func(x, result.x), label=u"拟合数据")
plt.legend(loc="lower center")
用basinhopping()
全局优化函数找出正弦波的三个参数。
前两个参数和其他求最小值的函数一样:目标函数和初始值。
由于它是全局优化函数,因此初始值的选择不太重要。
niter
是全局优化算法的的代次数,迭代的次数越多,就越有可能找到全域最优解
minimizer_kwargs
参数决定了所采用的局与最小值算法以及传递给此函数的参数。