python最小二乘法拟合实例

1）曲线拟合的最小二乘法
⑴两种逼近概念：
①插值：在节点处函数值相同；
②拟合：在数值上误差最小；
⑵最小二乘法主要解决的是插值中出现的三个问题：
① 高多项式会龙格现象；
② 定函数关系不合理(数据误差)；
⑶评估逼近效果：
假设“逼近”规律的近似函数为y=f(x), 即有yi﹡= f(xi)(i=1,2,…,m).它与观测值yi之差 δi = yi﹡ - yi = f(xi) – yi，i = 1, 2, …,m。这称为残差。残差大大小可以作为衡量近似函数好坏的标准。按照使残差的平方和∑δi²最小的规则求得近似函数y=f(x)的方法称为最佳平方逼近，也称为曲线拟合的最小二乘法。
举例：
1）假定某天的气温变化记录如下，使用最小二乘法找出这一天的气温变化规律.
t/h 0 1 2 3 4 5 6 7 8 9 10 11 12
T/℃ 15 14 14 14 14 15 16 18 20 22 23 25 28
t/h 13 14 15 16 17 18 19 20 21 22 23 24
T/℃ 31 32 31 29 27 25 24 22 20 18 17 16
要求：编程实现考虑下列类型的拟合函数，计算误差平方和，并作图比较效果.
（1）二次、三次、四次多项式拟合函数;
（2）形如
计算各个多项式的二范数的平方、残差的方差和残差绝对值的均值如表1-1所示.

残差的方差的意义在于评估各个偏差的离散程度，避免只是部分拟合的效果比较好，而部分偏差又非常大。残差绝对值的均值的意义在于评估总体水平上残差的大小。通过比较我们可以发现总体上是四次多项式的拟合效果比较好。

接下我们分别作出各个多项式在样本数据上的拟合效果，结果如图1-1所示。

通过图1-1我们也可以看出四次项的拟合效果比较好。
但通过增加多项式的最高次数之后，测试发现随着次数增加拟合的效果越好，出现了过拟合的现象，所以目前的评价体系还有问题。参考机器学习中的常见防止过拟合的方法是加入L1或者L2正则化项。
《目前还在学习阶段有些问题还请谅解》

代码如下：

import matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties
import numpy as np
from scipy.optimize import leastsq

x = np.arange(0,25,1,dtype='float')
y = np.array([15,14,14,14,14,15,16,18,20,22,23,25,
     28,31,32,31,29,27,25,24,22,20,18,17,16], dtype='float')

# param:最小化的目标| *_fun:拟合函数
# 二次
param0 = [0, 0, 0]
def quadratic_fun(a, x):
    k1, k2, b = a
    return k1 * x **2 + k2 * x + b
# 三次
param1 = [0, 0, 0, 0]
def cubic_fun(a, x):
    k1, k2, k3, b = a
    return k1 * x ** 3 + k2 * x**2 + k3 * x + b
# 四次
param2 = [0, 0, 0, 0, 0]
def fpower_fun(a, x):
    k1, k2, k3, k4, b = a
    return k1 * x ** 4 + k2 * x**3 + k3 * x**2 + k4 * x + b
# 自定义函数
param3 = [0, 0, 0]
def myfuns(s, x):
    a, b, c = s
    return a * np.exp(-b*(x-c))
param4 = [0, 0, 0, 0, 0, 0, 0, 0]
def testfun(a, x):
    k1, k2, k3, k4, k5, k6, k7, b = a
    return  k7 * x ** 7 + k6 * x ** 6 + k1 * x ** 5 + k2 * x**4 + k3 * x**3 + k4 * x ** 2 + k5 * x + b
# 偏差
def dist(a, fun, x, y):
    return fun(a, x) - y
# 作图
myfont = FontProperties(fname="C:\Windows\Fonts\msyh.ttc")
plt.figure()
plt.title(u'某天的气温变化', fontproperties=myfont)
plt.xlabel(u't/h')
plt.ylabel(u'T/摄氏度', fontproperties=myfont)
# 坐标轴的范围xmin, xmax, ymin, ymax
plt.axis([0, 24, 10, 35])
plt.grid(True)
plt.plot(x, y, 'k.')


funs = [quadratic_fun, cubic_fun, fpower_fun, myfuns, testfun]
params = [param0, param1, param2, param3, param4]
colors = ['blue', 'red', 'black', 'green', 'yellow']
fun_name = ['quadratic_fun', 'cubic_fun', '4power_fun', 'a*exp(-b*(t-c))', 'testfun']

for i, (func, param, color, name) in enumerate(zip(funs, params, colors, fun_name)):
    var = leastsq(dist, param, args=(func, x, y))
    plt.plot(x, func(var[0], x), color)
    print('[%s] 二范数: %.4f, abs(bias): %.4f, bias-std: %.4f' % (name,
          ((y-func(var[0], x))**2).sum(),
          (abs(y-func(var[0], x))).mean(), 
          (y-func(var[0], x)).std())
    )
plt.legend(['sample data', 'quadratic_fun', 'cubic_fun', '4power_fun', 'a*exp(-b*(t-c))', 'fivpower_fun'], loc='upper left') 
plt.show()

参考书目：
[M] 郑继明，朱伟，刘勇等 数值分析