馬永剛,張啟敏,劉俊梅

(1.寧夏大學數(shù)學統(tǒng)計學院,寧夏銀川 750021)

(2.榆林學院數(shù)學與統(tǒng)計學院,陜西榆林 719000)

1　引言

隨機微分方程是近幾年熱門的數(shù)學學科,它是微分方程、動力系統(tǒng)及隨機分析相互交叉的學科.在實際應用中,隨機微分方程的參數(shù)一般是未知的,因此參數(shù)估計問題成為一個重要的研究方向.關于此類方程參數(shù)估計問題已有很多結(jié)果,見文獻[1–3,9]等.在文獻[1]中,Bishwal詳細介紹了參數(shù)估計的理論方法和技巧,并運用到各種隨機模型中.在文獻[3]中,作者基于采樣的時間序列討論了非線性隨機微分方程的參數(shù)估計問題.

1926年意大利數(shù)學家Volterra提出著名的Lotka-Volterra微分方程模型,該模型引起眾多生物學家、數(shù)學家極大的興趣.Lotka-Volterra模型所表現(xiàn)的生態(tài)學現(xiàn)象,在自然界中處處可見.例如,病蟲害的周期爆發(fā),農(nóng)作物的豐收年與低產(chǎn)年的周期循環(huán),動物之間的捕食與被捕食等,這些現(xiàn)象都是Lotka-Volterra模型變化規(guī)律的具體展示.Lotka-Volterra模型的一般形式為[4]

其中αi,βi,γi(i=1,2)均為常數(shù),α1與α2分別表示兩種群的內(nèi)稟增長率,β1與γ2分別表示種內(nèi)作用系數(shù),γ1與β2分別表示種間作用系數(shù).就其生態(tài)意義可分為三類:競爭模型(αi＞ 0,βi＜ 0,γi＜ 0,i=1,2)、互惠模型 (αi＞ 0,βi＜ 0,γi＞ 0,i=1,2)、捕食 -被捕食模型(αi＞ 0,βi＜ 0,γ1＜ 0,γ2＞ 0,i=1,2).基于此模型,各種形式的Lotka-Volterra生態(tài)模型被提出,其中一類是受外界隨機因素影響的Lotka-Volterra生態(tài)模型.下面研究白噪聲擾動的Lotka-Volterra競爭模型[5]

其中r1,r2,a11,a12,a21,a22均為正常數(shù),ω1(t),ω2(t)為白噪聲且相互獨立.許多學者對此模型進行了研究,見文獻[4,5,6]等.其中最重要的結(jié)果之一,種群隨機持久生存,需要滿足下面條件

其中

當上述條件變?yōu)?/p>

兩種群中至少有一個種群隨機滅絕.

本文主要研究隨機種群Lotka-Volterra競爭模型的參數(shù)估計,即對方程組(1.2)中參數(shù)r1,a11,a12,r2,a21,a22,σ1,σ2應用最小二乘法理論進行估計,得到較好的估計結(jié)果.主要結(jié)果是參數(shù)的點估計值和區(qū)間估計,同時獲得影響估計區(qū)間長度的主要因素.最后給出了數(shù)值模擬,結(jié)果表明了該方法的可行性與有效性.

2　回歸模型

其中方程組(2.1)可轉(zhuǎn)化為

其中

進一步方程組(2.2)表示為

方程組(2.3)是簡單的線性回歸模型,應用線性回歸理論的方法估計參數(shù),估計過程基于最小二乘法,使下面兩個式子的值越小越好

為計算方便,方程(2.3)可寫成分塊矩陣的形式

其中

在塊矩陣中,子矩陣分別為

3　點估計

由多元線性回歸理論的公式估計參數(shù)β和η:

記

故(3.1)式可進一步化簡為

為了獲得參數(shù)的置信區(qū)間,需要估計參數(shù)β和η的方差,方差估計公式為

估計得到,p為被估計參數(shù)個數(shù).由(3.1),(3.3)式可得

等式(3.4),(3.5)化簡得

定理3.1方程組(3.6)中分別為方程組(2.1)中參數(shù)的漸進無偏估計,即當n→∞,

證類似文獻[7]證明.利用分別估計(3.2)式中參數(shù)得到被估參數(shù)的方差

4　區(qū)間估計

如果σ2已知,由最小二乘法回歸理論[8],參數(shù)估計值?β,?η的各分量都滿足正態(tài)分布.當觀察值的數(shù)量n足夠大時,可由?σ2代替σ2,因此(1?α)的置信區(qū)間(CIs)分別為

下面給出具體數(shù)據(jù),進行參數(shù)的估計.

例1在模型(1.2)式中,我們選擇滿足種群持久生存的數(shù)據(jù)進行模擬,即兩種群隨機持久存在的情形,其它狀態(tài)的數(shù)據(jù)也可獲得類似結(jié)論.假設參數(shù)真實值分別為r1=1,a11=0.3,a12=0.2,σ1=0.04,r2=1.2,a21=0.3,a22=0.4,σ2=0.05, 初始值為x1(0)=1.5,x2(0)=2.應用這些參數(shù)的值,通過Euler-Maruyama方法離散,我們獲得三組模擬數(shù)值xk1(t),xk2(t),數(shù)量分別為(A)n=5000;(B)n=10000;(C)n=50000;保存這些數(shù)據(jù)為模擬數(shù)據(jù)集.對每一組數(shù)據(jù),給步長?t=0.02,0.05,0.1進行模擬,給出一些模擬結(jié)果,以便比較真實值與估計值.在表1–4中,樣本的大小以符號“Size n”表示,且在表中第一列給出.其中表1,2樣本大小從5000增加到50000,在每一組模擬數(shù)據(jù)中,分別以三種不同步長進行模擬.模擬數(shù)據(jù)表明參數(shù)的估計值與真實值沒有明顯區(qū)別,絕對誤差與樣本大小有關,與步長?t無關,且隨著樣本數(shù)量的增大,絕對誤差越來越小.

表1:(r1,a11,a12)的估計模擬結(jié)果(真實值:r1=1,a11=0.3,a12=0.2)

表2:(r2,a21,a22)的估計模擬結(jié)果(真實值:r2=1,a21=0.3,a22=0.2)

其次給出各個參數(shù)置信水平為0.95的置信區(qū)間的一些模擬結(jié)果.表3,4中,給出參數(shù)的估計區(qū)間及區(qū)間長度,模擬數(shù)據(jù)表明隨著樣本量的增大,置信水平長度越來越小.

表3:置信水平為0.95的參數(shù)r1,a11,a12置信區(qū)間模擬結(jié)果(真實值:r1=1,a11=0.3,a12=0.2)

表4:置信水平為0.95的參數(shù)r2,a21,a22置信區(qū)間模擬結(jié)果(真實值:r2=1,a21=0.3,a22=0.2)

5　結(jié)論

本文對隨機兩種群Lotka-Volterra競爭模型應用最小二乘法理論進行參數(shù)估計,獲得參數(shù)r1,a11,a12,r2,a21,a22,σ1,σ2的估計值及估計區(qū)間.通過觀察結(jié)果,得到影響置信區(qū)間長度的主要因素.結(jié)果表明,參數(shù)估計區(qū)間的長度隨著樣本數(shù)量n的增加而減小,不依賴步長?t的長度.例1給出了具體的數(shù)值模擬.極大似然估計與貝葉斯估計方法是另外兩種參數(shù)估計的方法,在將來的工作中,將應用這些方法到隨機兩種群Lotka-Volterra模型.

[1]Bishwal JPN.Parameter estimation in stochastic differential equations[M].Berlin:Springer,2008.

[2]Kristensen NR,Madsen H,Young PC.Parameter estimation in stochastic grey-box model[J].Automatica,2004,40(2):225–237.

[3]Timmer J.Parameter estimation in nonlinear stochastic differential equatuons[J].Chaos Sol.Fract.,2000,11(15):2571–2578.

[4]Zhu C,Yin G.On competitive Lotka-Volterra model in random environments[J].J.Math.Anal.Appl.,2009,357(1):154–170.

[5]Jiang D Q,Ji C Y,Li X,O’Regan D.Analysis of autonomous Lotka-Volterra competition systems with random perturbation[J].J.Math.Anal.Appl.,2012,390:582–595.

[6]王佳,丁潔麗.Logistic回歸模型中參數(shù)極大似然估計的二次下界算法及其應用[J].數(shù)學雜志,2015,35(6):1521–1532

[7]Pan J F,Alison G,David G,Mao X R.Parameter estimation for the stochastic SIS epidemic model[J].Stat.Infer.Stoch Proc.,2014,17(1):75–98.

[8]Rawlings JO.Applied regression analysis:a research tool[M].Belmont,CA:Wadsworth,1988.

[9]洪志敏,閆在在.Volterra積分方程的蒙特卡羅數(shù)值求解方法[J].數(shù)學雜志,2016,36(2):425–436.