bTodo在线任务管理

简洁易用的在线任务管理系统 ...

↑我也要推荐

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

发布时间:2011-06-20 04:19:35, 关注:+22152, 赞美:+428, 不爽:+13

本文标签: numpy scipy

原始出处: 补完计划

前段时间帮@littlemorning做论文,要编程实现经济学的模型,其中主要用最小二乘拟合来估算一些函数的参数。科学计算的活一般来说都会用matlab,不过那样庞大的东西不是我所喜欢的。于是乎转向Python,一直对Python挺感兴趣,但是米有机会写点有用的东西,这次正好借这机会,体会一下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还有非常多的功能,足以应付日常的科学计算需求,以后有时间再深入挖掘一下。

如果你觉得本站对你有帮助,欢迎向本站赞助 :P

使用支付宝捐赠

Copyright© Python4cn(news, jobs) simple-is-better.com, 技术驱动:powered by web.py 空间主机:Webfaction

版权申明:文章转载已注明出处,如有疑问请来信咨询。本站为 python 语言推广公益网站,与 python 官方没有任何关系。

联系/投搞/留言: en.simple.is.better@gmail.com 向本站捐赠