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

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

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

1、理论推导

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的。

 2、仿真模型

想利用最近正在尝试使用的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/τ）因此正规方程的解也是只有一个元素的矩阵。将其转换为标量后取倒获得最优估计。

3、结果评估

加入噪声后的信号以下图所示，与一般状况的实测波形吻合度很高。

 

图2 模拟产生的带有噪声的阶跃响应 

 对这些加入噪声的信号进行线性校订后获得进行最小二乘估计的信号yt为下图所示。其基本趋势是一条斜率为（1/τ）的直线，和我预计的结果同样。

 

图3 对图2进行线性校订后的待估计信号

 

但能够明显的看到，因为（3）式左侧在V(t)的大小接近Vcc时对加入的白噪声进行了放大。所以当t递增时，由白噪声形成的信号的不肯定性大大增长了。也就是在套用正规方程时，t值较大时的噪声对结果的影响大于t值较小时的噪声对结果的影响。这是使用非线性校订函数须要付出的重要代价。

另外，屡次运行以上代码的获得 都是一个略小于实际τ（=0.1）的数值——约为0.099左右，也就是说， 不是无偏估计。这应该是因为线性校订函数（（3）式左侧），对于噪声noise大于0和小于0的部分进行了非对称的变换形成的。这虽然形成的误差不大，但也是使用非线性校订函数须要付出的代价。 

4、Jupyter notebook

上述练习的一个重要目的是练习使用Jupyter notebook，并在其中内嵌具备很好交互性的公式等信息。如下是部分程序运行效果的截图，虽然我对markdown语法还不熟悉，格式和排版还不够漂亮，但仍是可以明显的提升代码的可读性。

 

 

 图4 Jupyter notebook运行效果截图

 其中须要重点记录下的是公式代码的嵌入过程：

1）我首先在纸上手写了一些公式，用手机对其拍照，如： 

 

图5 手写的公式

  2）用mathpix tools对这些照片截图，并扫描（mathpix tools有windows版和iOS版，都可免费试用）。Mathpix能够直接获得LaTeX格式的公式表达式。下图是iOS版本的mathpix界面截图。

图6 iOS版本的mathpix截图

3）在Jupyter notebook上部的下拉菜单中选择单元格的格式为Markdown；将LaTeX格式的公式表达式粘贴到该单元格内，并在LaTeX公式表达式的先后加上“$$”符号，运行该单元格便可获得图4所示效果的公式了。

 

图7 在Jupyter notebook中输入LaTeX公式

 

 5、进一步的实际测试

工做作到这里离其实并无完，还应该作一个简单的实际电路实测一下。我会在后续的假期中抽半天时间，在STM32开发板上完成这个工做：用GPIO产生一个节约信号后，连续采集5个τ时间长度的信号，并代入正规方程求解，欢迎你们继续关注更新。

……