用 Python 做科学计算之最小二乘

这次我使用的是numpy和scipy这两个库，大家可以去官网了解一下,学习资料就是官方的文档了，写的还是非常详细的，英文不错的同学建议还是看看官方文档比较好。这两个库的优秀的中文资料还是挺少的，我只找到这一份比较好的文档，是介绍Python在科学计算领域内的应用的，说了很多内容，基本涵盖了科学计算的各个主要方面，其中提到numpy和scipy的地方说的比较简单，但是清晰，适合快速入门，作者貌似要出书，大家可以填一份问卷调查(文档主页上有链接)，可以得到最新的pdf文档。总的来说，numpy提供了一套适合数值计算的数据结构，而scipy提供了类似matlab库风格的函数。值得一提的是scipy的核心部分主要是用Fortran编写的，速度和可靠性都有很大保证。

好了，回到最小二乘的问题上，最小二乘拟合属于优化问题，在scipy的optimize子函数库中，提供了leastsq函数用于实现最小二乘。我们来看一个最简单的例子，假设需要拟合的函数是y = ax + b, 给定的一组数据,我们首先定义待拟合函数，然后定义一个函数求出真实值和拟合值之间的误差，然后调用leastsq函数进行拟合，请看代码：

import numpy as np
from scipy.optimize import leastsq

#待拟合的函数，x是变量，p是参数
def fun(x, p):
    a, b = p
    return a*x + b

#计算真实数据和拟合数据之间的误差，p是待拟合的参数，x和y分别是对应的真实数据
def residuals(p, x, y):
    return fun(x, p) - y

#一组真实数据，在a=2, b=1的情况下得出
x1 = np.array([1, 2, 3, 4, 5, 6], dtype=float)
y1 = np.array([3, 5, 7, 9, 11, 13], dtype=float)

#调用拟合函数，第一个参数是需要拟合的差值函数，第二个是拟合初始值，第三个是传入函数的其他参数
r = leastsq(residuals, [1, 1], args=(x1, y1))

#打印结果，r[0]存储的是拟合的结果，r[1]、r[2]代表其他信息
print r[0]

运行之后，拟合结果是

[2. 1.]

但是在这次实际的使用过程中，我拟合的函数不是这样简单的，其中的一个难点是待拟合函数是一个分段函数，需要判断自变量的值，然后给出不同的函数方程式，举个例子, 这样一个分段函数:当x > 3时，y = ax + b, 当x <= 3 时，y = ax – b, 用Python代码写一下：

def fun(x, p):
    a, b = p
    if (x > 3):
        return a*x + b
    else:
        return a*x - b

如果我们还是使用原来的差值函数进行拟合，会得到这样的错误：

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

原因很简单，我们现在的fun函数只能计算单个值了，如果传入的还是一个array，自然就会报错。那么怎么办呢？我也很郁闷，于是在scipy的maillist里寻求帮助, 外国牛牛们都很热心，很快就指出了问题。其实是我对于差值函数理解错了，leastsq函数所要传入的差值函数需要返回的其实是一个array, 于是我们可以这样修改差值函数：

def residuals(p, x, y):
    temp = np.array([0,0,0,0,0,0],dtype=float)
    for i in range(0, len(x)):
        temp[i] = fun(x[i], p)
    return temp - y

这样修改之后，我们就可以正常拟合啦。

其实scipy还有非常多的功能，足以应付日常的科学计算需求，以后有时间再深入挖掘一下。