1.线性函数拟合
先来看最简单的例子,对于线性函数
y=ax+b,可以使用最小二乘法来确定a、b两个系数。对于N个点,x值对应函数值与y值的差值平方和(残差平方和)作为目标函数,分别对a、b求偏导,令偏导为0,得到两个方程,联立解得系数。
设有N个点坐标为
(xi,yi),i=1,2,...,N,则残差平方和为
∑i=1N(axi+b−yi)2,展开来就是
(ax1+b−y1)2+(ax2+b−y2)2+...+(axN+b−yN)2,分别对a、b求偏导并令偏导为0得两个方程:
-
∑i=1N(axi+b−yi)xi=0
-
∑i=1N(axi+b−yi)=0
联立以上两个方程即可求解。
2.可化成线性函数的曲线拟合
幂函数拟合、双曲线(在一三象限的那种形式)拟合、指数函数拟合、对数函数拟合,都可以利用最小二乘法来解决。但是这些函数的偏导数可能很复杂,列出的方程组是非线性方程组,解算较为困难。因此我们需要通过变换点坐标把这些曲线转化为线性函数。
幂函数的一般方程为
y=axb,两边取对数可变换为
lny=blnx+lna,换元可得
y′=bx′+a′。对每个点(x,y)的x值和y值取对数可得到(x’,y’),化为线性函数拟合问题。
双曲线的一般方程为
y=ax+b1,两边取倒数可变换为
y1=ax+b,换元得
y′=ax+b。对每个点(x,y)的y值取倒数即可得到(x,y’),化为线性函数拟合问题。
指数函数的一般方程为
y=aebx,两边取自然对数可变换为
lny=lna+bx,换元得
y′=a′+bx。对每个点(x,y)的y值取对数即可得到(x,y’),化为线性函数拟合问题。
对数函数的一般方程为
y=alnx+b,换元得
y=ax′+b。对每个点(x,y)的x值取对数即可得到(x’,y),化为线性函数拟合问题。
3.系数间为线性关系的曲线拟合
再来看个例子
y=aex+bsinx+cx3,虽然这不是一个线性函数,也不能通过坐标变换化成线性函数,但如果设x为常量,a、b、c为变量,则a、b、c三个系数之间是线性关系的,因此残差平方和求偏导后得到的方程组依然是线性方程组。
设有N个点坐标为
(xi,yi),i=1,2,...,N,则残差平方和为
∑i=1N(aexi+bsinxi+cxi3−yi)2,展开来就是
(aex1+bsinx1+cx13−y1)2+(aex2+bsinx2+cx23−y2)2+...+(aexN+bsinxN+cxN3−yN)2,分别对a、b、c求偏导并令偏导为0得三个方程:
-
∑i=1N(aexi+bsinxi+cxi3−yi)exi=0
-
∑i=1N(aexi+bsinxi+cxi3−yi)sinxi=0
-
∑i=1N(aexi+bsinxi+cxi3−yi)xi3=0
联立以上三个方程即可求解。
典例:n阶多项式拟合
n阶多项式拟合同样采用最小二乘法。大体思路是对于N个点,用n+1个待定系数设出函数,将残差平方和作为目标函数,分别对n+1个系数求偏导,令偏导为0,得到n+1个方程,联立解得系数。
设有N个点坐标为
(xi,yi),i=1,2,...,N,设n次函数
f(x)=∑k=0nakxk,则其残差平方和为
∑i=1N∑k=0n(akxik−yi)2,展开来就是
(a0+a1x1+...+anx1n−y1)2+(a0+a1x2+...+anx2n−y2)2+...+(a0+a1xN+...+anxNn−yN)2,设其对aj求偏导并令偏导等于0,得
2(a0+a1x1+...+anx1n−y1)x1j+2(a0+a1x2+...+anx2n−y2)x2j+...+2(a0+a1xN+...+anxNn−yN)xNj=0,整理得
∑i=1N∑k=0nakxij+k−∑i=1Nyixij=0,j=0,1,...,n,可见j有n+1个取值,将n+1个点代入可以列出n+1个线性方程,利用线性代数知识可以解出系数值。
4.n阶多项式拟合实例
本实例演示使用python进行n阶多项式拟合,线性函数拟合其实就是一阶多项式拟合,而非线性函数拟合实例则放到下一期进行讲解。
1999-2004年长江排污量如下表:
年份 |
1995 |
1996 |
1997 |
1998 |
1999 |
2000 |
2001 |
2002 |
2003 |
2004 |
排污量/亿吨 |
174 |
179 |
183 |
189 |
207 |
234 |
220 |
256 |
270 |
285 |
用二次多项式来拟合以上数据,python代码如下:
import numpy as np
from matplotlib import pyplot
num=int(input())
x=[]
y=[]
for i in range(0,num):
x.append(int(input()))
y.append(int(input()))
z1 = np.polyfit(x, y, 2)
z = np.polyval(z1, x)
print(z1)
pyplot.plot(x,y,'.')
pyplot.plot(x,z,'-')
pyplot.show()
得到模型为
y=0.84x2−3349.94x+3336464.85,图像如下: