最小二乘法曲线拟合及Python实现
[本文PDF下载]
在各种信号处理中，如果某种检测到的信号和输入信号之间可以建立起一种相对稳定的对应关系，那我们就可以依据这种输入输出之间，构建一种转换关系；如果更进一步可用函数的方式表达出来，则似乎离真理就更靠近了一些。尤其对于一些工程应用，有这些对应关系则意味着可以投入应用。
图-1 模拟的输入输出采集信号示意图
图-1示意图中，如果横坐标是代表某种输入信号，纵坐标表示对应的输出信号。这如果要当作传感器，有无可能投入应用？有难度，看情况。主要是同一个输出信号，可能有两个以上的输入信号与之对应，意味着这个器件在所见的信号范围内输入和输出信号之间存在不确定性。但是如果将图示的输入信号缩小范围到其中的某一段斜率较大的单调区间，则可能成为一款很好的检测应用设备。
曲线拟合存在不同的方法。如果使用最小二乘法曲线拟合，如果暂时脱离了纯粹应用目的，对于上面的示意图会表示毫无压力。本文主要针对最小二乘法的基本公式推导，简化以及程序模拟展开。但在很多基本原理上，将不做太多笔墨。
-
最小二乘法的基本处理步骤
针对以上的数据，即存在以下的对应的输入和输出表[1]：
| 输入 |  |
 |
 |
 |
… |  |
 |
 |
 |
| 输出 |  |
 |
 |
 |
… |  |
 |
 |
 |
由于输入输出的未知对应关系，因此需要假定x和y之间存在以下的关系：
(1)
事实上，即使这样设定，除非理想情况，否则针对每个输入变量x_i{i∈(0,1,…n)}，其结果不会每个等于表[1]中的y_i。如果这个方程恰好与输入输出吻合，会有以下的结果。
我们简化为：
(2)
然而这种方程明显看出[X]是n×(r+1)的矩阵，离直接解算[A]还缺少一些条件，甚至是无解。所以转而求其次，针对当前每个输入，将拟合公式的输出y和实际测量值y_i值之间存在偏差尽量减小就算达成。最小二乘法的要求就是让这个拟合函数给出的结果y和相应的测量值y_i的偏差平方和最小时，则认为在采用公式（1）中的拟合公式达到最佳。所以就有以下公式：
(3)
式（3）我们也可以理解为两个n维向量 和 在空间上的距离达到最小。将式（1）带入式（3）得：
(4)
在式（4）中，
和
是已知量，反而作为原先的系数a_r是待求的未知数。如果让式（4）多元函数有最小值，则其对必要条件就是R对其中的每个变量的偏导等于0，其中
)。
（5）
整理“化简”该方程有：
（6）
到这里，我们似乎把问题搞得复杂化了。有必要将式（6）分解一下，方便我们理解。
如果将式（6）简述为：
可以看到，这里的[(1+r) x (1+r)]的方阵[Q]，其实就是：
即：
而方程右侧的[Z]，其实就是：
所以我们得到式（6）就是在式（2）两侧同时左乘矩阵
的转置矩阵
：
(7)
总之，我们要求算的拟合系数就在矩阵[A]中，根据式（7），可得：
(8)
殊途同归啊。有了以上的方程，只要
不是奇异矩阵，方程就有唯一解。实际分析处理方法中，式（7）反而是一种最直接的解决方式。之所以将方程（6）分解为式（7），主要是为了后续使用Python中已有的对矩阵操作的功能函数，否则对着式（6）中的一个个乘积和，无从下手。图-1中的拟合还等着结果。
-
Python实现最小二乘法曲线拟合
有了上面的分析，借助于Python已有的算法函数，处理最小二乘法的曲线拟合还是比较方便的。实际处理中，因为要提供模拟信号，输出图形，以及拟合公式等输出的功能，有一些额外的代码。这些功能都放到了一个类中。这个拟合功能提供了一个拟合阶数可指定的方式，默认阶数为1的线性拟合。
使用这个曲线拟合类的时候，一般分为以下几个步骤：
1. 生成二维的模拟数据，分别代表输入信号x轴和输出信号y轴
- 实际使用时，则可以将实测数据组传递给函数即可
- 在模拟生成数据的时候，借用随机数，给x，y两组数据都产生一定的噪声干扰
- X轴的输入信号在噪声干扰之后，如果要看到与之吻合的图形输出，则需要如类中的排列操作，但是排列与否和最终的拟合结果没有影响。
2. 生成拟合参数矩阵[A]
- 由输入矩阵[x]，生成中间的转换矩阵[X]
- 借用nump生成方阵
- 借用nump计算结果矩阵
- 借用nump计算参数矩阵[A]
3. 输出曲线拟合的结果方程（9阶拟合）
Function(x)=9.798298367046073+4.227514414846013x+17.260341567511098x**2-22.3475762650598x**3-128.91239180578455x**4-46.45306126416519x**5+152.35732377108258x**6+126.51471748420587x**7-53.97163898194201x**8-67.97621051238404x**9
4. 指定范围和分辨率输出曲线拟合后的图形（可以使用默认值）
图-2 拟合后输出曲线与原始信号示意图（左侧3阶拟合，右侧侧9阶拟合）
图-3 与理想输入输出信号比较
在安费诺传感器一些典型的传感器类型中，模拟输出的类型也比较多，尤其经典的诸如NTC，对于其信号的处理，拟合的操作无处不在。即使在我们数字调理的各类诸如压力，温湿度，振动等各种传感器，拟合也在其中担当着重要角色。
以上所述的代码可以在当前页面免费下载。
import matplotlib.pyplot as plt
import math
import numpy
import random
class LeastSqauresCF():
def __init__(self):
self.order = 1
self.Xvalues = None
self.Yvalues = None
self.fig = plt.figure()
self.subfig = self.fig.add_subplot(111)
#========================================
# Draw the input points/curve in a figure
#========================================
def ReadAnddrawInputPoints(self,Inputx=None,InputY=None, order=1):
try:
if Inputx is None or InputY is None:
return False
if(len(Inputx)!=len(InputY)):
return False
if order <=0:
return False
if type(order) !=int:
self.order = int(order)
else:
self.order = order
self.Xvalues = Inputx
self.Yvalues = InputY
self.subfig.plot(Inputx,InputY,color='m',linestyle='',marker='*')
return True
except Exception as e:
print(e)
return False
#========================================
# Calculate curve fitting factor matrix
#========================================
def produceFittingCurveFactors(self):
try:
#(1) Generate matrix [X]
matX=[]
for i in range(0,len(self.Xvalues)):
matx1=[]
for j in range(0,self.order+1):
dx=1.0
for l in range(0,j):
dx = dx * self.Xvalues[i]
matx1.append(dx)
matX.append(matx1)
#(2) Generate matrix [X]T.[X]
matX_Trans = numpy.matrix(matX).T
matX_FinalX = numpy.dot(numpy.matrix(matX_Trans),numpy.matrix(matX))
#(3) Generate matrix Y' =[X]T.[Y]
matFinalY = numpy.dot(matX_Trans,numpy.matrix(self.Yvalues).T)
#(4) Solve the function:[A] = [[X]T.[X]]**(-1).[X]T.[Y]
matAResult=numpy.linalg.solve(numpy.array(matX_FinalX),numpy.array(matFinalY))
return matAResult
except Exception as e:
print(e)
return None
#========================================
# Output fitting curve function
#========================================
def outputFittingCurveFunction(self, inputMatFactors=None):
if inputMatFactors is None:
return False
i = 0
strFitting="Function(x)="
for a in inputMatFactors:
#print(a[0])
if i==0:
strFitting +=str(a[0])
else:
strFitting +=("+"if a[0]>0 else"") +str(a[0])+(("x**"+str(i)) if i>1 else "x")
i+=1
print(strFitting)
return strFitting
#========================================
# draw the curve based on the result function
#========================================
def drawFittedCurve(self, xRangeMin=None,xRangeMax=None,matAResult=None,resolution=0.01):
try:
if matAResult is None:
return False
if xRangeMin is None or xRangeMax is None:
xRangeMin = self.Xvalues[0]
xRangeMax = self.Xvalues[-1]
#print('xRangeMin: ',xRangeMin, 'xRangeMax: ',xRangeMax)
xxa= numpy.arange(xRangeMin,xRangeMax,resolution)
yya=[]
for i in range(0,len(xxa)):
yy=0.0
for j in range(0,self.order+1):
dy=1.0
#x[i]**j
for k in range(0,j):
dy*=xxa[i]
#a[j]*(x[i]**j)
dy*=matAResult[j]
yy+=dy
yya.append(yy)
#print(xxa,yya)
self.subfig.plot(xxa,yya,color='g',linestyle='-',marker='')
self.subfig.legend()
plt.show()
return True
except Exception as e:
print(e)
return False
#===============================================
# main
#===============================================
if __name__=="__main__":
LS = LeastSqauresCF()
#==============================================
#(1-1) Generate simulation data set
#==============================================
x = numpy.arange(-1,1,0.02)
y = [(5*a+2)*(3*a*a+1)*numpy.sin(a*2)*numpy.cos(a*4)+5 for a in x]
#==============================================
#(1-2) Add some noises to input and output data
#==============================================
x_noised=[]
y_noised=[]
i=0
for yy in y:
xx = x[i]
Ns=float(random.randint(90,110))/100
x_noised.append(xx*Ns)
y_noised.append(yy*Ns+5*Ns)
i+=1
#=============================================
#(1-3) Sort x, y in case x's sequence was disturbed during adding noise
#=============================================
zip_x_y = zip(x_noised,y_noised)
sorted_zip = sorted(zip_x_y, key=lambda x:x[0])
sorted_x, sorted_y = zip(*sorted_zip)
#print(sorted_x,sorted_y)
#=============================================
#(2) Processing curve fitting
#=============================================
if LS.ReadAnddrawInputPoints(sorted_x,sorted_y,5):#plot the orginal data curve
MatrixFactor = LS.produceFittingCurveFactors()
LS.outputFittingCurveFunction(MatrixFactor)
LS.drawFittedCurve(matAResult=MatrixFactor)
