衛(wèi)澤林，趙 恒

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

0 引言

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

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

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

1 兩步估計(jì)基礎(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時的觀測樣本點(diǎn)為：

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

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

Rosenblatt于1955年提出：對每個點(diǎn)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)的估計(jì)。由此看出Rosenblatt方法與傳統(tǒng)直方圖方法的區(qū)別在于：分割區(qū)間隨著要估計(jì)的點(diǎn)x變化，x始終位于區(qū)間的中心位置，從而能夠獲得較好的效果。

引入函數(shù)：

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

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

如果K()

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

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

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

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

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

2 Lotka-Volterra模型的兩步估計(jì)

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

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

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

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

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

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

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

圖1 標(biāo)準(zhǔn)的x1和其估計(jì)值的x11的對比

圖2 標(biāo)準(zhǔn)的x2和其估計(jì)值的x22的對比

圖3 標(biāo)準(zhǔn)的x3和其估計(jì)值的x33的對比

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

3 核光滑時窗寬的選擇

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

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

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

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

本文使用Lotka-Volterra模型來進(jìn)行實(shí)驗(yàn)，模型中觀測時間為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àn)，求估計(jì)參數(shù)的均值，選擇該模型下的最優(yōu)窗寬。

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

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

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

4 用boostrap法估計(jì)兩步估計(jì)的方差

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

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

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

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

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

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

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

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

表3 boostrap法做兩步估計(jì)的方差估計(jì)

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

5 總結(jié)

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