加载中…
个人资料
渔婆的愚公
渔婆的愚公
  • 博客等级:
  • 博客积分:0
  • 博客访问:99,784
  • 关注人气:14
  • 获赠金笔:0支
  • 赠出金笔:0支
  • 荣誉徽章:
相关博文
推荐博文
谁看过这篇博文
加载中…
正文 字体大小:

Scipy,numpy,pylab包做最小二乘法拟合,并且作图

(2012-11-11 09:29:09)
标签:

python

scipy

numpy

pylab

it

分类: python

SciPy函数库在NumPy库的基础上增加了众多的数学、科学以及工程计算中常用的库函数。例如线性代数、常微分方程数值求解、信号处理、图像处理、稀疏矩阵等等。由于其涉及的领域众多、本书没有能力对其一一的进行介绍。作为入门介绍,让我们看看如何用SciPy进行插值处理、信号滤波以及用C语言加速计算。
假设有一组实验数据(x[i], y[i]),我们知道它们之间的函数关系:y = f(x),通过这些已知信息,需要确定函数中的一些参数项。例如,如果f是一个线型函数f(x) = k*x+b,那么参数k和b就是我们需要确定的值。如果将这些参数用 p 表示的话,那么我们就是要找到一组 p 值使得如下公式中的S函数最小:
S(\mathbf{p}) = \sum_{i=1}^m [y_i - f(x_i, \mathbf{p}) ]^2

这种算法被称之为最小二乘拟合(Least-square fitting)。

scipy中的子函数库optimize已经提供了实现最小二乘拟合算法的函数leastsq。下面是用leastsq进行数据拟合的一个例子:

 


# -*- coding: utf-8 -*- 

import numpy as np 

from scipy.optimize import 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) 

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进行数据拟合 

# 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()


这个例子中我们要拟合的函数是一个正弦波函数,它有三个参数 Aktheta ,分别对应振幅、频率、相角。假设我们的实验数据是一组包含噪声的数据 x, y1,其中y1是在真实数据y0的基础上加入噪声的到了。


通过leastsq函数对带噪声的实验数据x, y1进行数据拟合,可以找到x和真实数据y0之间的正弦关系的三个参数: A, k, theta。下面是程序的输出:



>>> 真实参数: [10, 0.34000000000000002, 0.52359877559829882] 

>>> 拟合参数 [-9.84152775  0.33829767 -2.68899335]

_images/scipy_intro_01.png

我们看到拟合参数虽然和真实参数完全不同,但是由于正弦函数具有周期性,实际上拟合参数得到的函数和真实参数对应的函数是一致的。

0

阅读 评论 收藏 转载 喜欢 打印举报/Report
  • 评论加载中,请稍候...
发评论

    发评论

    以上网友发言只代表其个人观点,不代表新浪网的观点或立场。

      

    新浪BLOG意见反馈留言板 电话:4000520066 提示音后按1键(按当地市话标准计费) 欢迎批评指正

    新浪简介 | About Sina | 广告服务 | 联系我们 | 招聘信息 | 网站律师 | SINA English | 会员注册 | 产品答疑

    新浪公司 版权所有