常微分方程兩步估計的方差估計

衛澤林，趙 恒

（西安電子科技大學a.數學與統計學院；b.通信工程學院，西安 710126）

0 引言

常微分方程理論在很多領域都有廣泛的應用，通過分析系統所建立的常微分方程模型，可以對客觀過程提供一個很好的描述。常微分方程的應用性通常依賴于某些參數，了解這些參數對研究常微分方程刻畫的動力系統是至關重要的，然而這些參數在實踐中往往是未知的，它們必須從帶有噪聲的觀測數據中推斷或估計出來[1-3]。

最小二乘估計具備良好的特性，因此在實際研究中有很大的作用，但是當方程為一個非線性的高維系統的時候，最小二乘估計就會有很大的局限性。這種情況下，需要在非線性高維的參數空間下找到使得函數達到最小的值，其解決過程通常是通過梯度的算法完成的。雖然這些算法的改進對參數估計有了很好的結果，但是在實際的應用中，仍需要很大的計算量才可以完成，這樣就會影響在實際應用中的效率。為了避免以上兩種算法的困難，提出了兩步估計法[4，5]，該方法的第一步只需要用核光滑估計[6]找出相對應的非參數估計，因此很大程度上節省了計算時間；第二步最小化方程殘差平方和就可以得到參數的估計值。

同時，參數的方差估計是一個很重要的指標。如果方差估計不可知，則難以進行很多其他的統計推斷，比如假設檢驗和區間估計等。但是針對兩步估計的方差估計沒有理論的結果。本文通過對一個三維的Lotka-Volterra種間競爭模型[7]用兩步估計進行參數估計，針對估計出來的參數用boostrap方法做方差估計。

1 兩步估計基礎知識

考慮下面的常微分方程系統[8]：

其中X(t)=(X1(t)，…，Xm(t))T∈Rm是m維的狀態變量，θ=(θ1，…，θd)T∈Θ，Θ∈Rd表示的是d維的未知參數，X?(t)是X(t)的導函數。

假設對第i個狀態變量Xi，在時間tj時的觀測樣本點為：

其中εij是服從獨立分布的模型隨機誤差，且滿足Eεij=0和 Var(ε)=σ2。

第一步：X(t)和X?(t)的估計可以通過非線性方法，例如核光滑，樣條函數[9]等實現，本文中利用核光滑[10]來估計。

Rosenblatt于1955年提出：對每個點x，Δx表示以x為中心，長為h的區間，即[x-h/2，x+h/2] 。以x1，...，xn落入 Δx的個數除以nh，即{j:xj∈Δx，j=1，...，n}/nh作為x密度函數f(x)的估計。由此看出Rosenblatt方法與傳統直方圖方法的區別在于：分割區間隨著要估計的點x變化，x始終位于區間的中心位置，從而能夠獲得較好的效果。

引入函數：

則上述f(x)的估計可表示為：

從Rosenblatt估計看出，由于W(x)是均勻分布，為估計f在點x的取值f(x)，與x距離h2以外的樣本點不起作用，而在距離h2以內的樣本點作用相同?？紤]到用整個樣本x1，...，xn估計f(x)，直觀上，越接近x的樣本點對f(x)的估計影響應當越大。Parzen于1962年把Rosenblatt估計推廣為核光滑估計：

如果K()

μ是一個給定的概率密度，則f(x)的核估計為：

其中K(μ)稱為核函數，h為窗寬。

可以根據以上核估計函數估計出X(t)和X?(t)，記光滑后的為X(t)和X?(t)分別為a(t)和b(t)。

第二步：估計參數θ的值。通過解下面的優化問題就可以得到參數θ的估計值：

其中H(θ)=(b(t)-F(a(t)，θ))2dt。

2 Lotka-Volterra模型的兩步估計

Logistic方程始于上世紀上半葉做昆蟲養殖實驗提出的“蟲口方程”，同一時代即在人口研究、魚群捕撈研究乃至生態系統研究上迅速發展起了Logistic方程的理論與應用研究，至今方興未艾。Logistic方程起源于生物領域的研究，但該方程可以說是提出了描述一切客觀系統發展過程的一個最簡單、最基本、最廣泛且最有用的數學模型。由于一維的Logistic方程僅限于描述單個事物或系統的發展過程，為了描述種群間的競爭關系Lotka和Volterra早在20世紀20年代，對一維Logistic方程進行直接擴展，提出了種群間競爭模型[11]。

考慮以下三維的Lotka-Volterra（LV）種間競爭模型：

將光滑過后的X(t)和X?(t)分別記為：

未知參數：

以第一個方程為例，記參數θ1=(β11，β12，β13)T的系數為h1(t)=(X1(1-X1/K)，X1X2，X1X3)，取K=1，則求解以下優化問題就可以得到θ1的估計：

對上式求導并令其為0，可解得最優解為：

類似的可以估計出其余兩個方程中的參數。

在做實驗時，利用模擬數據，設定參數為：

θ=(0 .01，0.01，0.01，0.02，0.02，-0.02，0.01，0.02，0.03)T，該值可以作為參數評價的一個標準。利用模擬數據做100次實驗，試驗中三個變量的噪聲均滿足正態分布N(0.0.022)。

表1 Lotka-Volterra模型中參數的估計結果

從表1中整體來看，各個參數的估計值都比較準確，多次實驗后求得的參數估計的均值也比較準確，ARE較小。圖1至圖3中，x1、x2和x3分別是狀態變量的標準值，x11、x22和x33分別是將估計出來的參數代入方程，再用龍格庫塔解出的狀態變量的值，兩者相比較，從圖中可以看出，兩步估計效果很好。

圖1 標準的x1和其估計值的x11的對比

圖2 標準的x2和其估計值的x22的對比

圖3 標準的x3和其估計值的x33的對比

但是在第一步估計中，用光滑的核函數估計X(t)時，其中需要選擇合適的窗寬，不同的窗寬對參數估計的結果影響較大，因此表1的估計結果不一定是最優的，下文將介紹該模型中最優窗寬的選擇。

3 核光滑時窗寬的選擇

在核光滑估計問題中，窗寬的選擇[12]直接影響到估計的好壞，而核函數形式對估計的影響甚微。窗寬h越大，則光滑越好,但可能忽視有用信息，而擬合不好。反之，h越小，則擬合較好，但可能光滑不夠，無法把有用信息和干擾分開。

窗寬參數的選擇有許多種方法，通常第一步的核光滑估計可以采用數據驅動的交叉驗證或者建立在漸近平均整合平方誤差基礎上的替換法，其最優窗寬為：

Liang和Wu[13]在2008年提出另一種兩步估計的最優窗寬為：

其中an是一個收斂于0但收斂速度慢于log-1n的序列，在本文他們選擇an=log-1/16n。

本文使用Lotka-Volterra模型來進行實驗，模型中觀測時間為1～20天，則n=20，由于 log20＞1，因此an＜1，那么在這個模型中有：

本文分別取an1=log-1/106n，an2=log-1/16n，an3=log-1/6n，對應窗寬分別為hn1=0.4244，hn2=0.3967，hn3=0.3538，每個窗寬下做100次實驗，求估計參數的均值，選擇該模型下的最優窗寬。

表2 Lotka-Volterra模型中不同窗寬下估計參數的ARE

由表2可以看出，無論哪個參數，其變化規律和結果都是一樣的：

其結果可以看出，對于Lotka-Volterra模型，窗寬的選擇使用hopt=n-1/5時，估計的參數更加準確，而若使用公式hn=n-2/7an,當觀測值個數一定時，窗寬會受到限制，只能取到一定小范圍內的值，而不一定可以取到最優窗寬，從以上模型的實驗來看，使用hn=n-2/7an公式時，在可取得的窗寬范圍內，窗寬越大，ARE越小，估計越準確。

4 用boostrap法估計兩步估計的方差

本文用Lotka-Volterra模型做實驗時都用的模擬數據，因此可以模擬多組實驗數據以計算方差，但是在實際問題中通常只會有一組數據，也只能根據這一組數據估計參數，難以計算方差。Boostrap是現代統計學較為流行的一種統計方法，在小樣本時效果很好[14，15]。通過方差的估計可以構造置信區間等，其運用范圍得到進一步延伸。所以本文使用boostrap估計兩步估計參數的方差。

這種方法的過程比較簡潔：

（1）首先采用重抽樣技術從原始樣本中抽取一定數量的樣本，這個過程可以重復。而在兩步估計中，目標是解決:

這是一個積分問題，可以從積分區間入手，將積分區間[a，b]均勻的等分為N份，看成N個原始樣本，從中可放回的抽取N個樣本區間。

（2）根據抽出的樣本區間計算常微分方程中的參數的估計值：

（3）重復上述過程n次，得到參數的n個估計值。

（4）最后計算上述參數n個估計值的樣本方差，這樣就得到了常微分方程估計參數的方差估計。

下面本文將這種方法運用在Lotka-Volterra模型中，用20個觀測值做實驗，試驗中第一步估計的窗寬選擇第3部分中得到的最優窗寬值hopt=n-1/5。用上述方法計算可以估計出來參數的方差。而由于實驗中是模擬數據，因此也可以計算出其標準方差，由于方差數值較小，比較不方便，而用標準差比較性質也不會發生改變，因此表3中使用標準差比較實驗結果。

表3 boostrap法做兩步估計的方差估計

表3中，第二列是由多次實驗計算得到兩步估計的參數估計值，再計算這些參數估計值標準的方差，進而求出標準差，第三列是由用boostrap方法取樣，估計方差并求出標準差?？梢钥闯?，用boostrap做常微分方程參數方差估計的方法計算出的方差與標準方差相比相對誤差較小，這種方法做方差估計較精確可行。

5 總結

常微分方程模型在很多領域中都有廣泛的應用，比如生物學和醫學等。而常微分方程中參數的估計方法也多種多樣。本文利用基于光滑的核函數兩步估計法來估計參數，這種方法比傳統的參數估計方法更簡潔有效，計算效率高。其次本文研究了兩步估計中窗寬對估計參數方差的影響。最后基于最優窗寬，提出了用boostrap法做兩步估計的方差估計。本文通過對一個三維的Lotka-Volterra種間競爭模型進行參數估計，針對估計出來的參數做方差估計，并驗證了該方法的有效性。方差的估計對很多統計推斷問題很有意義，例如假設檢驗，區間估計等，本文研究了兩步估計的方差估計，但是由于兩步估計本身的一些性質，其中很多因素都會對估計結果造成影響，因此很難將所有影響因素考慮進去，進而找到一個最優的方差估計。這些需要在未來的工作中繼續研究。其次本文提到的方法在高維常微分方程上的推廣也有待進一步研究。