Python基础教程学习第九日：SciPy-数值计算库

本文参考：http://old.sebug.net/paper/books/scipydoc/scipy_intro.html

最小二乘法

最小二乘法：

# -*- coding: utf-8 -*-
import numpy as np
from scipy.optimize import leastsq  #最小二乘法拟合函数leastsq
import pylab as 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)   # 以0为起点，-2*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进行数据拟合，p为需要拟合的参数
# residuals为计算误差的函数
# p0为拟合参数的初始值
# args为需要拟合的实验数据
plsq = leastsq(residuals, p0, args=(y1, x))

print u"真实参数:", [A, k, theta]
print u"拟合参数", plsq[0] # 实验数据拟合后的参数

pl.plot(x, y0, label=u"真实数据")
pl.plot(x, y1, label=u"带噪声的实验数据")
pl.plot(x, func(x, plsq[0]), label=u"拟合数据")
pl.legend()
pl.show()

p0里放的是A、k、theta的初始值，这个值可以随意指定。往后随着迭代次数增加，A、k、theta将会不断变化，使得residuals函数的值越来越小。

func函数里指出了待拟合函数的函数形状。

residuals函数为误差函数，我们的目标就是不断调整A、k、theta使得residuals不断减小。这里的residuals函数和神经网络中常说的cost函数实际上是一回事，只不过这里更简单些而已。

必须注意一点，传入leastsq函数的参数可以有多个，但必须把参数的初始值p0和其它参数分开放。其它参数应打包到args中。

leastsq的返回值是一个tuple，它里面有两个元素，第一个元素是A、k、theta的求解结果，第二个元素我暂时也不知道是什么意思，先留下来。

输出：

>>> 真实参数: [10, 0.34000000000000002, 0.52359877559829882]
>>> 拟合参数 [ 9.57495101  0.34221818  0.55335328]

函数最小值

# -*- coding: utf-8 -*-
# 本程序用各种fmin函数求卷积的逆运算

import scipy.optimize as opt
import numpy as np

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

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

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, h_等数列，调用各种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
h_ = np.random.rand(n)

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

if __name__ == "__main__":
test_n(200, 20, 0.1)

convolve_func(h)为需要求最小值的函数，h为该函数的参数，也是需要求解最小值的参数。

fminfunc(convolve_func, h_)， fminfunc为函数变量，在这里为：opt.fmin、opt.fmin_powell、opt.fmin_cg、opt.fmin_bfgs，也就是分别调用这些最小化函数。第一个变量为需要计算最小值的函数，第二个变量为该函数的变量的初始值。

输出为：

fmin
ーーーーーーーーーーー
error of y: 0.00568756699607
error of h: 0.354083287918

fmin_powell
ーーーーーーーーーーー
error of y: 0.000116114709857
error of h: 0.000258897894009

fmin_cg
ーーーーーーーーーーー
error of y: 0.000111220299615
error of h: 0.000211404733439

fmin_bfgs
ーーーーーーーーーーー
error of y: 0.000111220251551
error of h: 0.000211405138529

非线性方程求解

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
print f(result)

fsolve实际上求解的是f=0这个函数，即：5*x1 + 3 = 0；4*x0*x0 - 2*sin(x1*x2) = 0；x1*x2 - 1.5 = 0。[1,1,1]为初始值。

fsolve函数在调用函数f时，传递的参数为数组，因此如果直接使用数组中的元素计算的话，计算速度将会有所降低，因此这里先用float函数将数组中的元素转换为Python中的标准浮点数，然后调用标准math库中的函数进行运算。

在对方程组进行求解时，fsolve会自动计算方程组的雅可比矩阵，如果方程组中的未知数很多，而与每个方程有关的未知数较少时，即雅可比矩阵比较稀疏时，传递一个计算雅可比矩阵的函数将能大幅度提高运算速度。雅各比矩阵如下：


使用雅可比矩阵的fsolve实例如下，计算雅可比矩阵的函数j通过fprime参数传递给fsolve，函数j和函数f一样，有一个未知数的解矢量参数x，函数j计算非线性方程组在矢量x点上的雅可比矩阵。由于这个例子中未知数很少，因此程序计算雅可比矩阵并不能带来计算速度的提升，如下：

# -*- coding: utf-8 -*-
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
]

def j(x):
x0 = float(x[0])
x1 = float(x[1])
x2 = float(x[2])
return [
[0, 5, 0],
[8*x0, -2*x2*cos(x1*x2), -2*x1*cos(x1*x2)],
[0, x2, x1]
]  # 这里参考雅各比矩阵的计算方式，为各个函数对各个变量的偏导数

result = fsolve(f, [1,1,1], fprime=j)
print result
print f(result)

输出：

[-0.70622057 -0.6        -2.5       ]
[0.0, -9.1260332624187868e-14, 5.3290705182007514e-15]

插值函数

interpolate库提供了许多对数据进行插值运算的函数。

# -*- coding: utf-8 -*-
import numpy as np
import pylab as pl
from scipy import interpolate

x = np.linspace(0, 2*np.pi+np.pi/4, 10)
y = np.sin(x)

x_new = np.linspace(0, 2*np.pi+np.pi/4, 100)
f_linear = interpolate.interp1d(x, y)
tck = interpolate.splrep(x, y)
y_bspline = interpolate.splev(x_new, tck)

pl.plot(x, y, "o",  label=u"原始数据")
pl.plot(x_new, f_linear(x_new), label=u"线性插值")
pl.plot(x_new, y_bspline, label=u"B-spline插值")
pl.legend()
pl.show()

除了以上interp1d、splrep之外，还有其他的插值方式，详见：https://my.oschina.net/u/2306127/blog/600628