最小二乘法(least squares)的曲线拟合(curve fitting)_最小二乘法拟合曲线

第三十八篇 最小二乘法的曲线拟合

如果我们想得到一个通过大量由实验或者计算机程序获得的数据点的函数,它实际是在寻找一个“最适合”数据的函数,而不是一个完全经过所有点。可以采用各种策略来最小化各个数据点之间的误差和逼近函数。其中最著名的是最小二乘法，它让用户能够自由选择在曲线拟合中使用的函数类型。这种方法也被称为具有两个或两个以上的自变量的“线性回归”或“多元线性回归”，。

最小二乘法

考虑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
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)
    
    
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

sparin
1

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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

spabac
1

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]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

终端输出结果如下