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

衛(wèi)澤林，趙 恒

（西安電子科技大學a.數(shù)學與統(tǒng)計學院；b.通信工程學院，西安 710126）

0 引言

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

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

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

1 兩步估計基礎(chǔ)知識

考慮下面的常微分方程系統(tǒng)[8]：

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

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

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

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

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

引入函數(shù)：

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

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

如果K()

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

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

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

第二步：估計參數(shù)θ的值。通過解下面的優(yōu)化問題就可以得到參數(shù)θ的估計值：

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

2 Lotka-Volterra模型的兩步估計

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

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

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

未知參數(shù)：

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

對上式求導(dǎo)并令其為0，可解得最優(yōu)解為：

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

在做實驗時，利用模擬數(shù)據(jù)，設(shè)定參數(shù)為：

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

表1 Lotka-Volterra模型中參數(shù)的估計結(jié)果

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

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

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

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

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

3 核光滑時窗寬的選擇

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

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

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

其中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，對應(yīng)窗寬分別為hn1=0.4244，hn2=0.3967，hn3=0.3538，每個窗寬下做100次實驗，求估計參數(shù)的均值，選擇該模型下的最優(yōu)窗寬。

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

由表2可以看出，無論哪個參數(shù)，其變化規(guī)律和結(jié)果都是一樣的：

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

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

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

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

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

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

（2）根據(jù)抽出的樣本區(qū)間計算常微分方程中的參數(shù)的估計值：

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

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

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

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

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

5 總結(jié)

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