非线性函数的最小二乘拟合——兼论Jupyter notebook中使用公式 [原创]

Posted 2021-03-22 helesheng

tags:

篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了非线性函数的最小二乘拟合——兼论Jupyter notebook中使用公式 [原创]相关的知识，希望对你有一定的参考价值。

突然有个想法，利用机器学习的基本方法——线性回归方法，来学习一阶RC电路的阶跃响应，从而得到RC电路的结构特征——时间常数τ（即R*C）。回答无疑是肯定的，但问题是怎样通过最小二乘法、正规方程，以更多的采样点数来降低信号采集噪声对τ估计值的影响。另外，由于最近在捣鼓Jupyter和numpy这些东西，正好尝试不用matlab而用Jupyter试试看。结果是意外的好用，尤其是在Jupyter脚本中插入LaTeX格式的公式的功能，真是太方便了！尝试了直接把纸上手写的公式转换到Jupyter脚本中的常见工具软件。

以下原创内容欢迎网友转载，但请注明出处： https://www.cnblogs.com/helesheng

一、理论推导

1．线性回归分析及正规方程

传统意义说，线性回归问题是用最小二乘法（即正规方程），解决线性方程组的均方误差最小化问题。已知输出输入X是由多个变量构成的行向量，W是系数向量（列向量），b为偏置

在机器学习中，把每次的输入x作为一行组成更大的矩阵，即每一行代表一个样本，该矩阵称为设计矩阵X^(train)。若样本数为k，则X^(train)有k行，每个样本都会得到一个输出y，将y集合成一个列向量Y^(train)，k个相同的b也组成列向量b。为简化表达，将b简化为附加在系数列向量W最后的常数b，X^(train)则在每行的最后增加一个1，W（列向量）的最后增加一个待估计的b。为了使估计的结果：

和Y^(train)之差的平方和最小，有正规方程可以求解W:

技术图片

2.一阶RC电路的阶跃响应

一阶RC电路的电路模型如下图所示。

技术图片

图1 一阶RC电路

幅度为Vcc的阶跃信号从Vin处输入，在Vout处测量输出。解微分方程可得自变量为时间t的响应函数。

技术图片

其中时间常数τ = R*C。我希望通过测量阶跃信号输入条件下，实际RC电路的响应曲线V(t)，并通过V(t)来估计时间常数τ。如果做纯理论推导，只要知道任意时刻t0的输出电压V(t0)就可以解方程（2）得到τ。但在实际电路中电压V(t0)的测量往往受到诸多干扰的影响，并不准确。是否可以测量多个t值时刻对应的V(t)，并利用正规方程（1）得到一个统计意义上最优的估计呢？是接下来要解决的问题。

3.非线性函数的最小二乘估计

仔细观察适用正规方程的目标函数（0）式的特点，可以发想让非线性的要让（2）式能够使用正规方程，必须让：

1）含有待估计的变量τ的函数充当（0）式中的系数W，而设计矩阵X^(train)则可以由含有时间t或测量电压V(t)的函数充当。

2） W和X^(train)之间的关系必须是简单的相乘。

显然，只有用时间t的序列充当设计矩阵X^(train)，才有可能使W和X^(train)之间的关系必须是相乘。至于其他的非线性部分则可以通过等效变换，变换到等式的另一侧来充当Y^(train)。综上，可以将（2）式变换为（3）式。

技术图片

（3）式的整个左边可以计算得到Y^(train)；（3）式右边的时间t的序列在构成设计矩阵X^(train)，1/τ则构成相当于（0）式中的系数矩阵W。这样就可以通过正规方程（2）式来求解统计意义上最优的估计了。当然，从（3）式也可以看出，经过线性校正的模型中系数向量W只有一个元素——是个标量，偏置b也是恒等于0的。

二、仿真模型

想利用最近正在尝试使用的Jupyter和numpy替代熟悉的Matlab，验证上述非线性函数最小二乘估计的想法。下面先建立一个模型：

1）输入为幅度Vcc为3.3V的阶跃信号；

2）时间常数τ为0.1秒；

3）简单模拟采样间隔的随机性：是间隔等于delta1=0.0015秒和delta1=0.0011秒的两个序列的叠加。

4）采样总长度为n=5倍τ；

5）信号上叠加了幅度约为20mV的白噪声——至于为什么是20mV，将在后续部分详细介绍。

利用python和numpy进行数值仿真的代码如下：

 1 import numpy as np
 2 import matplotlib.pyplot as plt
 3 tao=0.1#时间常数
 4 vcc=3.3#电源电压
 5 n=5#时长取时间常数tao的n倍
 6 delta1=0.0015#第一种采样间隔
 7 delta2=0.0011#第一种采样间隔
 8 t1=delta1*np.arange(np.ceil(n*tao/delta1))
 9 t2=delta2*np.arange(np.ceil(n*tao/delta2))
10 t=np.append(t1,t2)#为演示最小二乘拟合的功能，故意结合两种采样率下的时间点
11 t.sort()#对t进行排序
12 plt.plot(t)
13 s=vcc*(1-np.exp(-t/tao))#理论的波形曲线
14 plt.plot(t,s)#注意这里的plot函数使用了x轴和y轴两个轴，因为s中的数据不是均匀的
15 N_amp=np.exp(-n)*vcc
16 N_amp
17 noise = np.random.uniform(-N_amp, N_amp, (len(t)))#噪声：正负np.exp(-5)*3.3之间均匀分布
18 s_nr=s+noise#加入噪声后的信号
19 plt.plot(t,s_nr)
20 yt=np.log(vcc/(vcc-s_nr))
21 plt.plot(t,yt)
22 yt=np.mat(yt)
23 yt=yt.T
24 x=np.mat(t)#将X转换为矩阵数据类型
25 x=x.T#正规方程中的x应该是个列向量
26 w=(np.linalg.inv(x.T*x))*x.T*yt#求解正规方程
27 E_tao = np.round(1/float(w),4)#对时间常数的tao的估计，保留到4位小数
28 E_tao

非线性函数的最小二乘拟合

说明：

1）其间序列包含了两个等差数列t1和t2的融合，它们的间隔互质，没有重复，目的是模拟采样时间的随机性。最后用sort()方法又对时间序列进行排序的目的是为了后续容易观察波形更直观。如果仅仅为了使用正规方程，是不需要重新排序的。

2）校正的非线性方程（3）和原始方程（2）相比有一个重大的缺陷就是：左侧求对数的括号内的数值不能为负——如果是纯理论推导，这是不可能发生的。但假如随机噪声后的V(t)是有可能大于阶跃幅度Vcc的，此时括号内将变为一个负数，使得计算变得没有意义。我的解决之道是将假如的随机噪声幅度限制在仿真截止时刻V(t)和Vcc之差的范围内，及代码中的N_amp。由于仿真的结束时刻为n（=5）个τ，所以N_amp的值等于np.exp(-n)*vcc。
这样做没有任何理论依据，仅仅是受限于（3）和（2）式之间的非完全等价变换——属于线性化校正需要付出的代价。

3）由于待估计的参数只有一个（1/τ）所以正规方程的解也是只有一个元素的矩阵。将其转换为标量后取倒得到最优估计。