最小二乘法(least squares)的曲线拟合(curve fitting) |
您所在的位置:网站首页 › 最小二乘法所有公式 › 最小二乘法(least squares)的曲线拟合(curve fitting) |
第三十八篇 最小二乘法的曲线拟合
如果我们想得到一个通过大量由实验或者计算机程序获得的数据点的函数,它实际是在寻找一个“最适合”数据的函数,而不是一个完全经过所有点。可以采用各种策略来最小化各个数据点之间的误差和逼近函数。其中最著名的是最小二乘法,它让用户能够自由选择在曲线拟合中使用的函数类型。这种方法也被称为具有两个或两个以上的自变量的“线性回归”或“多元线性回归”,。 最小二乘法考虑np数据点(x1, y1),(x2, y2),…,(np, ynp),其中,x表示自变量[x1, x2,…],xnv]T, y为因变量。自变量的数目nv可以取任何期望的值,虽然在传统的线性回归中它被设为1。所需的“最好匹配”的函数可以写成这样的形式 使用最小二乘法去推导最好的拟合直线,对于np个数据点(x1, y1),(x2, y2),…,(xnp , ynp )。 这个例子只有一个自变量x,因此nv=1。而且,如果一个线性方程能拟合数据,下面的方程牵扯到两个未知常量(k=2) 为了使用最小二乘法,所提出的曲线拟合函数必须是方程(5.39)的“线性”标准形式。在工程分析中经常遇到的一些有用的曲线拟合函数最初并不是这种形式,但是可以通过简单的变量变换转化为标准形式。 例如,考虑尝试采取“幂数定律”形式的函数 一个阻尼振子有一个固有频率w=91.7/s。在自由振动中,测量出随着时间振幅衰减的数据为 分为一个主程序,和两个子程序,分别为天际线存储向量的乔列斯基分解子程序sparin,和逆向迭代求解的子程序spabac。详情可见以天际线存储矩阵系数的乔列斯基分解 #最小二乘法的曲线拟合 import numpy as np import math import B npt=5;nv=1;nc=2 kdiag=np.zeros((nc),dtype=np.int64) kv=np.zeros(nc*(nc+1)//2) f=np.zeros((nc)) c=np.zeros((nc)) x=np.zeros((npt,nv)) y=np.zeros((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(c,np.transpose(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]终端输出结果如下 |
今日新闻 |
点击排行 |
|
推荐新闻 |
图片新闻 |
|
专题文章 |
CopyRight 2018-2019 实验室设备网 版权所有 win10的实时保护怎么永久关闭 |