非線性函數的最小二乘擬合——兼論Jupyter notebook中使用公式 [原創]

時間 2019-12-20

標籤非線性函數最小二乘擬合 jupyter notebook 使用公式原創简体版

原文原文鏈接

忽然有個想法，利用機器學習的基本方法——線性迴歸，來學習一階RC電路的階躍響應，從而獲得RC電路的結構特徵——時間常數τ（即R*C）。回答無疑是確定的，但問題是怎樣經過最小二乘法、正規方程，以更多的採樣點數來下降信號採集噪聲對τ估計值的影響。另外，因爲最近在搗鼓Jupyter和numpy這些東西，正好嘗試不用matlab而用Jupyter試試看。結果是意外的好用，尤爲是在Jupyter腳本中插入LaTeX格式的公式的功能，真是太方便了！嘗試了直接把紙上手寫的公式轉換到Jupyter腳本中的常見工具軟件。
python

如下原創內容歡迎網友轉載，但請註明出處： https://www.cnblogs.com/heleshengwindows

1、理論推導markdown

1．線性迴歸分析及正規方程app

傳統意義說，線性迴歸問題是用最小二乘法（即正規方程），解決線性方程組的均方偏差最小化問題。已知輸出輸入X是由多個變量構成的行向量，W是係數向量（列向量），b爲偏置dom

在機器學習中，把每次的輸入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:ide

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運行效果截圖

其中須要重點記錄下的是公式代碼的嵌入過程：