SciPy —— （一）func & optimize

最新推荐文章于 2020-12-08 08:44:26 发布

Artoria____

最新推荐文章于 2020-12-08 08:44:26 发布

阅读量352

点赞数

分类专栏： NLP # SciPy

本文链接：https://blog.csdn.net/artoria_qzh/article/details/102153206

版权

NLP 同时被 2 个专栏收录

13 篇文章 0 订阅

订阅专栏

SciPy

1 篇文章 0 订阅

订阅专栏

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$ 是概率统计中常出现的一个特殊函数：

$\Gamma(z)=\int_{0}^{\infty}t^{z-1}e^{-t}dt$

在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))

在这里插入图片描述
该过程实现了对方程组： $\begin{bmatrix} 5x+3=0\\ 4x_0^2-2sin(x_1x_2)=0\\ x_1x_2-1.5=0 \end{bmatrix}$ 的求解

在对方程组进行求解时，fsolve()会自动计算方程组在某点对各个未知数变量的偏导数，这些偏导数组成一个二维数组，也就是雅可比矩阵。
如果方程组中的未知数很多，而与每个方程有关联的未知数较少，即雅可比矩阵比较稀疏时，将计算雅可比矩阵的函数作为参数传给fsolve()，能大幅度提高计算速度。

（二）最小二乘拟合

假如有一组数据(xi, yi)且假设它们满足某种函数关系yi = f(xi)。

若f()=kx+b，为了确定k和b的值，我们可以用最小二乘拟合来确定

如果用p表示函数中需要确定的参数，则目标是找到一组p使得函数s的值最小：
$S(P)=\sum_{i=1}^{m}[y_i-f(x_i,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)对x和y的偏导函数为：

$\frac{\partial f}{\partial x}=-2+2x-400x(y-x^2)$

$\frac{\partial f}{\partial y}=200y-200x^2$

而Newton-CG算法需要计算海森矩阵，它是一个由自变量为向量的实值函数的二阶偏导数构成的方块矩阵

则对于该f(x, y)来说，海森矩阵为一个二阶矩阵：

$\begin{bmatrix} 2(600x^2-200y+1)&-400x \\ -400x& 200 \end{bmatrix}$

最后用各种最小值优化算法计算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参数决定了所采用的局与最小值算法以及传递给此函数的参数。

Artoria____

关注

0
点赞
踩
4

收藏

觉得还不错? 一键收藏
0
评论
SciPy —— （一）func & optimize

Scipy在NumPy的基础上增加了众多的数学计算、科学计算以及工程计算中常用的模块。接下来会对其慢慢介绍目录一、常数和函数（一）常数（二）函数二、optimize（一）非线性方程组求解（二）最小二乘拟合一、常数和函数（一）常数Scipy的constants模块包含了众多的物理常数from scipy import constants as Cprint(C.c) ...
复制链接

扫一扫