最⼩⼆乘法（leastsquares）的曲线拟合（curvefitting）
第三⼗⼋篇最⼩⼆乘法的曲线拟合
如果我们想得到⼀个通过⼤量由实验或者计算机程序获得的数据点的函数，它实际是在寻⼀个“最适合”数据的函数,⽽不是⼀个完全经过所有点。可以采⽤各种策略来最⼩化各个数据点之间的误差和逼近函数。其中最著名的是最⼩⼆乘法，它让⽤户能够⾃由选择在曲线拟合中使⽤的函数类型。这种⽅法也被称为具有两个或两个以上的⾃变量的“线性回归”或“多元线性回归”，。 最⼩⼆乘法
考虑np数据点(x1, y1)，(x2, y2)，…，(np, ynp)，其中，x表⽰⾃变量[x1, x2，…]，xnv]T, y为因变量。⾃变量的数⽬nv可以取任何期望的值，虽然在传统的线性回归中它被设为1。所需的“最好匹配”的函数可以写成这样的形式 其中fj(x)， j = 1,2，…，k为选函数的属于x, Cj, j = 1,2，…，k是常数，通过最⼩⼆乘法过程得到的最优解。线性回归中的“线性”⼀词仅指上⾯⽅程中模型对Cj常数的依赖性。fj(属于x)函数由⽤户来控制，如果需要可以是⾮线性的。
⽬标是使F(x)尽可能接近y，因此，考虑以下每个np数据点上这些量之间的差的平⽅和
上⾯⽅程给出的误差项可以通过E对每个常数Cj的偏导来最⼩化，并使结果等于零，从⽽
这样的k个线性联⽴⽅程的对称形式可以写成矩阵形式
使⽤线性求解中的⽅法求出C1, C2，…。最后，将优化后的Cj常数代回曲线拟合⽅程1，根据需要进⾏进⼀步插值。
计算实例:
使⽤最⼩⼆乘法去推导最好的拟合直线，对于np个数据点(x1, y1),(x2, y2),…,(xnp , ynp )。
这个例⼦只有⼀个⾃变量x，因此nv=1。⽽且，如果⼀个线性⽅程能拟合数据，下⾯的⽅程牵扯到两个未知常量(k=2)
其中f1(x)=1，f2(x)=x
从⽅程1可以得到
它能够求解得到
拟合直线
把c1和c2代⼊到⽅程F(x)中，经典线性回归⽅程能获得。
数据的相关联系数可以得到，
决定系数由r2得到，是数据线性依赖程度的度量。
数据的线性化
为了使⽤最⼩⼆乘法，所提出的曲线拟合函数必须是⽅程(5.39)的“线性”标准形式。在⼯程分析中经常遇到的⼀些有⽤的曲线拟合函数最初并不是这种形式，但是可以通过简单的变量变换转化为标准形式。
例如，考虑尝试采取“幂数定律”形式的函数
上⾯的函数不是⼀个标准的形式，因为常量B作为幂数，不是简单乘法的系数。通过取两边的对数得到
现在如果做⼀个替换，X=lnx和Y=lny得到
或者⼀个标准形式为
其中C1=lnA和C2=B，变量的改变得到了标准形式下lnx和lny的线性关系。最⼩⼆乘法将得到最合适的系数C1和C2，代⼊到原来的⽅程中得到A=e**c1和B=C2
这个进程可以命名叫做‘线性化’过程，能被⽤在⼤量初始不为标准形式的函数中。⼀些例⼦中，直转化⼀个变量就能实现线性化。⼀些⽅程的线性化展⽰在下表中。
计算实例
⼀个阻尼振⼦有⼀个固有频率w=91.7/s。在⾃由振动中，测量出随着时间振幅衰减的数据为
拟合的⽬标曲线为：
使⽤最⼩⼆乘法去发现这个函数，然后计算临界阻尼ζ.
这个指数衰减曲线拟合函数并不是标准形式，所以需要转化。初始⽅程为，
参考之前的转化表能得到线性形式
因此
其中Y = lny X = t，C1 = lny0，C2 = - 91.7ζ。根据最开始的最⼩⼆乘法拟合函数，建⽴联⽴⽅程的函数为f1(X)=1, f2(X) = X。
对于⼿算时，推荐采⽤表格表⽰，得到
从⽅程1得到，
得到解为
因为C2 =−18.210，那么就得出ζ =−18.210/−91.7 = 0.20。
程序如下
分为⼀个主程序，和两个⼦程序，分别为天际线存储向量的乔列斯基分解⼦程序sparin，和逆向迭代求解的⼦程序spabac。详情可见
#最⼩⼆乘法的曲线拟合
import numpy as np
import math
import B
npt=5;nv=1;nc=2
s((nc),dtype=np.int64)
s(nc*(nc+1)//2)
s((nc))
s((nc))
s((npt,nv))
s((npt))
x[:,0]=(29,50,74,103,118)
y[:]=(1.6,23.5,38.0,46.4,48.9)
print('最⼩⼆乘法的曲线拟合')
my=sum(y)/npt
def f54(x,f):
f[0]=1.0
f[1]=math.log(x[0])
return f54
for i in range(1,npt+1):
f54(x[i-1,:],f);ic=0
for j in range(1,nc+1):
c[j-1]=c[j-1]+f[j-1]*y[i-1]
for k in range(1,j+1):
ic=ic+1;kv[ic-1]=kv[ic-1]+f[j-1]*f[k-1] for i in range(1,nc+1):
kdiag[i-1]=i*(i+1)/2
B.sparin(kv,kdiag)
B.spabac(kv,c,kdiag)
print('最优函数系数')
for i in range(1,nc+1):
print('{:13.4e}'.format(c[i-1]),end='') print()
print('数据点和拟合点')
print(' (x(i),i=1,',nv,'),y,yi')
sd=0;es=0
for i in range(1,npt+1):
f54(x[i-1,:],f);yi=np.dot(anspose(f)) sd=sd+(y[i-1]-my)**2;es=es+(y[i-1]-yi)**2 for j in range(1,nv+1):
print('{:13.4e}'.format(x[i-1,j-1]),end='') print('{:13.4e}'.format(y[i-1]),end='') print('{:13.4e}'.format(yi))
r2=(sd-es)/sd
print('相关联系数的平⽅',r2)
sparin
def sparin(kv,kdiag):
#对称天际线矩阵的乔列斯基分解
n=kdiag.shape[0]
kv[0]=kv[0]**0.5
for i in range(2,n+1):
ki=kdiag[i-1]-i
l=kdiag[i-2]-ki+1
for j in range(int(l),i+1):
x=np.float64(kv[ki+j-1])
kj=kdiag[j-1]-j
if j!=1:
ll=kdiag[j-2]-kj+1
ll=max(l,ll)
if ll!=j:
m=j-1
for k in range(int(ll),m+1):
x=x-np.float64(kv[ki+k-1]*kv[kj+k-1]) kv[ki+j-1]=x/kv[kj+j-1]
kv[ki+i-1]=x**0.5
spabac
def spabac(kv,loads,kdiag):
#天际线矩阵的乔列斯基前后迭代
n=kdiag.shape[0]
loads[0]=loads[0]/kv[0]
for i in range(2,n+1):
ki=kdiag[i-1]-i
l=kdiag[i-2]-ki+1
x=loads[i-1]
if l!=i:
m=i-1
for j in range(int(l),m+1):
x=x-kv[ki+j-1]*loads[j-1]
loads[i-1]=x/kv[ki+i-1]
for it in range(2,n+1):
i=n+2-it
ki=kdiag[i-1]-i
x=loads[i-1]/kv[ki+i-1]
loads[i-1]=x
l=kdiag[i-2]-ki+1
if l!=i:
m=i-1
for k in range(int(l),int(m+1)):
loads[k-1]=loads[k-1]-x*kv[ki+k-1]
loads[0]=loads[0]/kv[0]
终端输出结果如下
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议