1. 曲线拟合的最小二乘法及python实现
首页地址https://kz16.top/na
本文来源https://kz16.top/na/leastSquare
最小二乘法一般提法: 对于给定的一组数据(xi,yi),i=0,1,⋯,m,要求在函数空间Φ=span{ϕ0,ϕ1,⋯,ϕn}找到一个函数y=S∗(x),使误差加权平方和
∣∣δ∣∣22=i=0∑mω(xi)[S∗(xi)−yi]2=S(x)∈ϕmini=0∑mω(xi)[S(xi)−yi]2.
这里ω(x)>0是权函数,这里的
S(x)=a0ϕ0+a1ϕ1+⋯+anϕn(x),(n<m).(1)
其中a0,a1,⋯,an为待求系数。
由(1)式可知,最小二乘法可以进一步转化为求有限维多元函数
I(a0,a1,⋯,an)=i=0∑mω(xi)[j=0∑najϕj(xi)−f(xi)]2(2)
的极小点(a0∗,a1∗,⋯,an∗)问题。这里f(xi)=yi。根据多元函数极值的必要条件,
∂ak∂I=2i=0∑mω(xi)[j=0∑najϕj(xi)−f(xi)]ϕk=0,k=0,1,⋯,n.
上式其实就是一个关于a=(a0,a1,⋯,an)T的线性方程组,矩阵形式为
Ga=d,(3)
其中
G=[gkj]k=0,1,⋯,nj=0,1,⋯,n,gkj=i=0∑mω(xi)ϕk(xi)ϕj(xi)
以及
d=[dk]k=0,1,⋯,n,dk=i=0∑mω(xi)ϕk(xi)f(xi)
因为基函数ϕ0,ϕ1,⋯,ϕn线性无关,因此矩阵G是一个满秩矩阵,从而线性方程组(3)存在唯一解。 根据gkj=gjk,因此矩阵G是对称矩阵。
我们当然最希望线性方程组(3)中的矩阵G是一个对角矩阵,也就是需要满足
gkj=i=0∑mω(xi)ϕk(xi)ϕj(xi)={0,Ak>0,j=k,j=k.(4)
我们称满足上面(4)式基函数族{ϕj}j=0,1,⋯,n为关于点集{xi}i=0,1,⋯,m带权{w(xi)}i=0,1,⋯,m的正交基函数族。
如果采用这种带权正交基函数,则线性方程组(3)的解为
ak∗=Akdk=∑i=0mω(xi)ϕk2(xi)∑i=0mω(xi)ϕk(xi)f(xi)
且平方误差为
∣∣δ∣∣22=∣∣f∣∣22−k=0∑kAk(ak∗)2.
已知一组实验数据如下表,求它的拟合曲线
xi |
1 |
2 |
3 |
4 |
5 |
f(xi) |
4 |
4.5 |
6 |
8 |
8.5 |
频数 |
2 |
1 |
3 |
1 |
1 |
很显然,这里的频数可以当成权,这组实验数据的点分布在一条直线附近,因此可以选择线性拟合,基函数取
ϕ0=1,ϕ1=x
import numpy as np
def phi(x):
return np.array([1.0+0*x, x])
x = np.array([1, 2, 3, 4, 5])
y = np.array([4, 4.5, 6, 8, 8.5])
w = np.array([2, 1, 3, 1, 1])
phixsw = phi(x)*np.sqrt(w)
leftG = phixsw.dot(phixsw.T)
rightd = phi(x).dot(y*w)
ans = np.linalg.solve(leftG, rightd)
print("Left matrix G = ")
print(leftG)
print("Right vector G = ")
print(rightd)
print("Solver ans = ")
print(ans)
结果
Left matrix G =
[[ 8. 22.]
[22. 74.]]
Right vector G =
[ 47. 145.5]
Solver ans =
[2.56481481 1.2037037 ]
因此拟合曲线为
f(x)=2.56481481+1.2037037x
注1: 用python或matlab运算为了提高效率,可以把问题改成矩阵向量的运算,从而可以直接使用内置函数
首页地址https://kz16.top/na
本文来源https://kz16.top/na/leastSquare
仅供参考,部分漏洞在所难免
欢迎转载,转载请指明来源,
请勿用于商业
作者:Zhou Bingzhen
2019.11.21