python最小二乘法拟合曲线

zhangqinghao

最小二乘拟合
当面对一堆数据时，常见的一种处理方法即是拟合这些数据，做出一条与原曲线近似重合的曲线并求出方程。

下面这一组数据来源于《数值分析》--李庆扬。

t/min        1        2        3        4        5        6        7        8        9        10        11        12        13        14        15        [折叠] 16
y/*10-3        4.00        6.40        8.00        8.80        9.22        9.50        9.70        9.86        10.00        10.20        10.32        10.42        10.50        10.55        10.58        10.60
可以看到这些点在坐标轴所近似的图形为双曲线或指数形式，即：

双曲线型 1/y = a + b/x，即y = x /(ax + b)
指数型 y = a * eb/x
下面分别用这两种模型来拟合，最后计算最大误差和均方误差，获得合理的拟合曲线。

#!/usr/bin/env python
# -*- coding: UTF-8 -*-
'''通過一組數據擬合線性方程組并繪制圖形。'''

import numpy as np
from scipy.optimize import leastsq
import pylab as pl


def hyperbola_residuals(p, y, x):
    '''計算双曲线模型的誤差并返回。'''
    return y - x / (p[0] * x + p[1])


def exponent_residuals(p, y, x):
    '''計算指數模型的誤差并返回。'''
    return y - p[0] * np.exp(p[1] / x)


def hyper_value(x, p):
    '''根據變量和參數計算双曲线函數的值。'''
    return  x / (p[0] * x + p[1])


def exp_value(x, p):
    '''根據變量和參數計算指数函數的值。'''
    return p[0] * np.exp(p[1] / x)


# 猜測的初始參數值
p0 = [1, 1]

# 數值
x = np.arange(1, 17, 1)
y = np.array([4.00, 6.40, 8.00, 8.80, 9.22, 9.50, 9.70, 9.86, 10.00, 10.20, 10.32, 10.42, 10.50, 10.55, 10.58, 10.60])

# 擬合，調用參數分別為誤差residuals, 猜測初始值, 需要擬合的實驗數據
hyper_plsq = leastsq(hyperbola_residuals, p0, args=(np.reciprocal(y), np.reciprocal(x)))
exp_plsq = leastsq(exponent_residuals, p0, args=(y, x))
# 分別輸出兩組結果
print('双曲线函数参数：', hyper_plsq[0])
print('指数函数参数:', exp_plsq[0])

# 分別計算兩個模擬的均方根誤差
hyper_err_ufunc = np.frompyfunc(lambda y, x: hyperbola_residuals(hyper_plsq[0], np.reciprocal(y), np.reciprocal(x)), 2, 1)
hyper_errors = hyper_err_ufunc(y, x)
hyper_rms = np.sqrt(np.sum(np.square(hyper_errors)))
exp_err_ufunc = np.frompyfunc(lambda y, x: exponent_residuals(exp_plsq[0], y, x), 2, 1)
exp_errors = exp_err_ufunc(y, x)
exp_rms = np.sqrt(np.sum(np.square(exp_errors)))
print('两个曲线模型的均方根误差:', hyper_rms, exp_rms)

# 繪制原曲線及擬合後的曲線
pl.plot(x, y, 'b^-', label='Origin Line')
pl.plot(x, hyper_value(x, hyper_plsq[0]), 'gv--', label='Hyperbola Fitting Line')
pl.plot(x, exp_value(x, exp_plsq[0]), 'r*', label='Exponent Fitting Line')
pl.axis([0, 18, 0, 18])
pl.legend()
# Save figure
pl.savefig('scipy01.png', dpi=96)
输出：

双曲线函数参数： [ 2.23437501  1.76562499]
指数函数参数: [ 11.35718526  -1.07275843]
两个曲线模型的均方根误差: 0.413291610267 0.331399405496
由最后的rms可以看出，指数函数模型更加合适。最后的图像：

protected.png

SixPy

转发别人的帖子，要贴源地址
python代码要保持格式，否则没法看，更没法用！

看你发了这么多贴，都是这样的，你赶快修复格式，以及 给出源地址。
否则就要当灌水处理了~