数学杂志  2018, Vol. 38 Issue (2): 367-374   PDF    
扩展功能
加入收藏夹
复制引文信息
加入引用管理器
Email Alert
RSS
本文作者相关文章
马永刚
张启敏
刘俊梅
具有随机扰动的Lotka-Volterra竞争模型的参数估计
马永刚1,2, 张启敏1, 刘俊梅2    
1. 宁夏大学数学统计学院, 宁夏 银川 750021;
2. 榆林学院数学与统计学院, 陕西 榆林 719000
摘要:本文研究了两种群随机Lotka-Volterra竞争模型的参数估计的问题.利用最小二乘法,获得了点估计及(1-α)置信区间估计,同时得到了影响置信区间长度的因素.最后给出数值模拟,结果表明该方法的可行性与有效性.
关键词随机Lotka-Volterra竞争模型    回归模型    最小二乘法    点估计    置信区间    
PARAMETER ESTIMATION FOR LOTKA-VOLTERRA COMPETITION MODEL WITH RANDOM PERTURBATIONS
MA Yong-gang1,2, ZHANG Qi-min1, LIU Jun-mei2    
1. School of Mathematics and Statistics, Ningxia University, Yinchuan 750021, China;
2. School of Mathematics and Statistics, Yulin University, Yulin 719000, China
Abstract: In this paper, the problem of parameter estimation for two species stochastic Lotka-Volterra competition model is studied. By using the method of least squares, we obtain the point estimate and the (1 -α) confidence interval, and we get the factors which influence the confidence interval length. Finally, the numerical simulation results show that the method is feasible and effective.
Key words: Lotka-Volterra competition model     regression model     least squares method     point estimators     interval estimation    
1 引言

随机微分方程是近几年热门的数学学科, 它是微分方程、动力系统及随机分析相互交叉的学科.在实际应用中, 随机微分方程的参数一般是未知的, 因此参数估计问题成为一个重要的研究方向.关于此类方程参数估计问题已有很多结果, 见文献[1-3, 9]等.在文献[1]中, Bishwal详细介绍了参数估计的理论方法和技巧, 并运用到各种随机模型中.在文献[3]中, 作者基于采样的时间序列讨论了非线性随机微分方程的参数估计问题.

1926年意大利数学家V.Volterra提出著名的Lotka-Volterra微分方程模型, 该模型引起众多生物学家、数学家极大的兴趣. Lotka-Volterra模型所表现的生态学现象, 在自然界中处处可见.例如,病虫害的周期爆发, 农作物的丰收年与低产年的周期循环, 动物之间的捕食与被捕食等, 这些现象都是Lotka-Volterra模型变化规律的具体展示. Lotka-Volterra模型的一般形式为[4]

$ \begin{equation} \left\{ \begin{aligned} \frac{dx_1(t)}{dt} &= x_1(t)[\alpha_1+\beta_1x_1(t)+\gamma_1x_2(t)], \\ \frac{dx_2(t)}{dt} &= x_2(t)[\alpha_2+\beta_2x_1(t)+\gamma_2x_2(t)], \end{aligned} \right. \end{equation} $ (1.1)

其中$\alpha_i, \beta_i, \gamma_i(i=1, 2)$均为常数, $\alpha_1$$\alpha_2$分别表示两种群的内禀增长率, $\beta_1$$\gamma_2$分别表示种内作用系数, $\gamma_1$$\beta_2$分别表示种间作用系数.就其生态意义可分为三类:竞争模型$(\alpha_i>0, \beta_i <0, \gamma_i<0, i=1, 2)$、互惠模型$(\alpha_i>0, \beta_i<0, \gamma_i>0, i=1, 2)$、捕食-被捕食模型$(\alpha_i>0, \beta_i<0, \gamma_1<0, \gamma_2>0, i=1, 2)$.基于此模型, 各种形式的Lotka-Volterra生态模型被提出, 其中一类是受外界随机因素影响的Lotka-Volterra生态模型.下面我们研究白噪声扰动的Lotka-Volterra竞争模型[5]

$ \begin{equation} \left\{ \begin{aligned} dx_1(t) &= x_1(t)[r_1-a_{11}x_1(t)-a_{12}x_2(t)]dt+\sigma_1x_1(t)d\omega_1(t), \\ dx_2(t) &= x_2(t)[r_2-a_{21}x_1(t)-a_{22}x_2(t)]dt+\sigma_2x_2(t)d\omega_2(t), \end{aligned} \right. \end{equation} $ (1.2)

其中$r_1, r_2, a_{11}, a_{12}, a_{21}, a_{22}$均为正常数, $\omega_1(t), \omega_2(t)$为白噪声且相互独立.许多学者对此模型进行了研究, 见文献[4, 5, 6]等.其中最重要的结果, 种群随机持久生存, 需要满足下面条件

$ a_{11}>a_{12}>0, a_{22}>a_{21}>0, \sigma_1>0, \sigma_2>0, r_1>\frac{a_{12}}{a_{22}}r_2, \\ r_2>\frac{a_{21}}{a_{11}}r_1, r_1-\frac{\sigma_1^2}{2}>\frac{a_{12}}{a_{22}}(r_2-\frac{\sigma_2^2}{2}), r_2-\frac{\sigma_2^2}{2}>\frac{a_{21}}{a_{11}}(r_1-\frac{\sigma_1^2}{2}), \\ \frac{a_{21}x_1^*}{2}\sigma_1^2+\frac{a_{12}x_2^*}{2}\sigma_2^2<\min\{a_{21}(a_{11}-a_{12})(x_1^*)^2, a_{12}(a_{22}-a_{21})(x_2^*)^2\}, $

其中

$ x_1^*=\frac{r_1a_{22}-r_2a_{12}}{a_{11}a_{22}-a_{12}a_{21}}, x_2^*=\frac{r_2a_{11}-r_1a_{21}}{a_{11}a_{22}-a_{12}a_{21}}. $

当上述条件变为

$ \begin{aligned} r_1<\frac{\sigma_1^2}{2} \quad \text{或} \quad r_2<\frac{\sigma_2^2}{2}, \end{aligned} $

两种群中至少有一个种群随机灭绝.

本文主要研究随机种群Lotka-Volterra竞争模型的参数估计, 即对方程组(1.2)中参数$r_1, a_{11}, a_{12}, r_2, a_{21}, a_{22}, \sigma_1, \sigma_2$应用最小二乘法理论进行估计, 得到较好的估计结果.主要结果是参数的点估计值和区间估计, 同时获得影响估计区间长度的主要因素.最后给出了数值模拟, 结果表明了该方法的可行性与有效性.

2 回归模型

$ \big\{ \binom {x_1^{(k)}} {x_2^{(k)}} \big\}_{k=0}^n$是模型(1.2)的观测值, 给定初始值$\binom {x_1^{(0)}} {x_2^{(0)}}$和步长$\Delta t.$考虑在小区间$[k\Delta t, (k+1)\Delta t]$上, 利用Euler-Maruyama方法, 得到离散化形式

$ \begin{equation} \left\{ \begin{aligned} x_1^{(k+1)}-x_1^{(k)} &= x_1^{(k)}[r_1-a_{11}x_1^{(k)}-a_{12}x_2^{(k)}]\Delta t+\sigma_1x_1^{(k)}\Delta\omega_1^{(k)}, \\ x_2^{(k+1)}-x_2^{(k)} &= x_2^{(k)}[r_2-a_{21}x_1^{(k)}-a_{22}x_2^{(k)}]\Delta t+\sigma_2x_2^{(k)}\Delta\omega_2^{(k)}, \end{aligned} \right. \end{equation} $ (2.1)

其中$ \Delta\omega_1^{(k)}=\omega_1^{(k+1)}-\omega_1^{(k)}, \Delta\omega_2^{(k)}=\omega_2^{(k+1)}-\omega_2^{(k)}.$方程组(2.1)可转化为

$ \begin{equation} \left\{ \begin{aligned} y_1^{(k)} &= r_1\sqrt{\Delta t}-a_{11}x_1^{(k)}\sqrt{\Delta t}-a_{12}x_2^{(k)}\sqrt{\Delta t}+ \frac{\sigma_1\Delta\omega_1^{(k)}}{\sqrt{\Delta t}}, \\ y_2^{(k)} &= r_2\sqrt{\Delta t}-a_{21}x_1^{(k)}\sqrt{\Delta t}-a_{22}x_2^{(k)}\sqrt{\Delta t}+ \frac{\sigma_2\Delta\omega_2^{(k)}}{\sqrt{\Delta t}}, \end{aligned} \right. \end{equation} $ (2.2)

其中

$ y_1^{(k)}=\frac{x_1^{(k+1)}-x_1^{(k)}}{x_1^{(k)}\sqrt{\Delta t}}, y_2^{(k)}=\frac{x_2^{(k+1)}-x_2^{(k)}}{x_2^{(k)}\sqrt{\Delta t}}, \\\varepsilon_1^{(k)}=\frac{\sigma_1\Delta\omega_1^{(k)}}{\sqrt{\Delta t}}\sim N(0, \sigma_1), \varepsilon_2^{(k)}=\frac{\sigma_2\Delta\omega_2^{(k)}}{\sqrt{\Delta t}}\sim N(0, \sigma_2). $

进一步方程组(2.2)表示为

$ \begin{equation} \left\{ \begin{aligned} y_1^{(k)} &= r_1\sqrt{\Delta t}-a_{11}x_1^{(k)}\sqrt{\Delta t}-a_{12}x_2^{(k)}\sqrt{\Delta t}+\varepsilon_1^{(k)}, \\ y_2^{(k)} &= r_2\sqrt{\Delta t}-a_{21}x_1^{(k)}\sqrt{\Delta t}-a_{22}x_2^{(k)}\sqrt{\Delta t}+ \varepsilon_2^{(k)}, \end{aligned} \right. \end{equation} $ (2.3)

方程组(2.3)都是简单的线性回归模型, 应用线性回归理论的方法估计参数, 估计过程基于最小二乘法, 使下面两个式子的值越小越好

$ \min \big \{ \sum\limits_{k=1}^n (y_1^{(k)}-r_1\sqrt{\Delta t}-a_{11}x_1^{(k)}\sqrt{\Delta t}-a_{12}x_2^{(k)}\sqrt{\Delta t}) \big \}, \\ \min \big \{ \sum\limits_{k=1}^n (y_2^{(k)}-r_2\sqrt{\Delta t}-a_{21}x_1^{(k)}\sqrt{\Delta t}-a_{22}x_2^{(k)}\sqrt{\Delta t}) \big\}. $

为计算方便, 方程(2.3)可写成分块矩阵的形式

$ \begin{eqnarray} Y=X\theta+\varepsilon. \end{eqnarray} $ (2.4)

其中

$ \begin{eqnarray*} Y = \begin{pmatrix} y_1\\ y_2\\ \end{pmatrix}, \quad X = \begin{pmatrix} A_{n\times3} & O_{n\times3}\\ O_{n\times3} & A_{n\times3}\\ \end{pmatrix}, \quad \theta = \begin{pmatrix} \beta\\ \eta\\ \end{pmatrix}, \quad \varepsilon = \begin{pmatrix} \varepsilon_1\\ \varepsilon_2\\ \end{pmatrix} . \end{eqnarray*} $

在块矩阵中, 子矩阵分别为

$ \begin{eqnarray*} &y_1 = \begin{pmatrix} y_1^{(1)}\\ y_1^{(2)}\\ \vdots \\ y_1^{(n)}\\ \end{pmatrix}, \quad y_2 = \begin{pmatrix} y_2^{(1)}\\ y_2^{(2)}\\ \vdots \\ y_2^{(n)}\\ \end{pmatrix}, \quad A_{n\times3}= \begin{pmatrix} \sqrt{\Delta t} &\, -x_1^{(1)}\sqrt{\Delta t} &\, -x_2^{(1)}\sqrt{\Delta t} \\ \sqrt{\Delta t} &\, -x_1^{(2)}\sqrt{\Delta t} &\, -x_2^{(2)}\sqrt{\Delta t} \\ \vdots &\, \vdots &\, \vdots \\ \sqrt{\Delta t} &\, -x_1^{(n)}\sqrt{\Delta t} &\, -x_2^{(n)}\sqrt{\Delta t}\\ \end{pmatrix}, \quad \\& O_{n\times3}= \begin{pmatrix} 0 &\, 0 &\, 0 \\ 0 &\, 0 &\, 0 \\ \vdots &\, \vdots &\, \vdots \\ 0 &\, 0 &\, 0 \\ \end{pmatrix}, \varepsilon_1 = \begin{pmatrix} \varepsilon_1^{(1)}\\ \varepsilon_1^{(2)}\\ \vdots \\ \varepsilon_1^{(n)}\\ \end{pmatrix}, \quad \varepsilon_2 = \begin{pmatrix} \varepsilon_2^{(1)}\\ \varepsilon_2^{(2)}\\ \vdots \\ \varepsilon_2^{(n)}\\ \end{pmatrix}, \quad \beta = \begin{pmatrix} r_1\\ a_{11}\\ a_{12}\\ \end{pmatrix}, \quad \eta = \begin{pmatrix} r_2\\ a_{21}\\ a_{22}\\ \end{pmatrix}. \end{eqnarray*} $
3 点估计

由多元线性回归理论的公式估计参数$\beta\text{和}\eta:$

$ \begin{equation} \label{eq:1} \begin{aligned} \hat{\theta} &= \begin{pmatrix} \hat{\beta}\\ \hat{\eta}\\ \end{pmatrix} =(X^{T}X)^{-1}(X^{T}Y)\\ &=\Big[\begin{pmatrix} A^{T} & O\\ O & A^{T}\\ \end{pmatrix} \begin{pmatrix} A & O\\ O & A\\ \end{pmatrix}\Big]^{-1} \begin{pmatrix} A^{T}y_1\\ A^{T}y_2\\ \end{pmatrix}\\ &= \begin{pmatrix} (A^{T}A)^{-1}(A^{T}y_1)\\ (A^{T}A)^{-1}(A^{T}y_2)\\ \end{pmatrix}. \end{aligned} \end{equation} $ (3.1)

$ \begin{eqnarray*} B= \begin{vmatrix} n\Delta t &-\Delta t\sum x_1^{(k)} &-\Delta t\sum x_2^{(k)}\\ -\Delta t\sum x_1^{(k)} &\Delta t\sum (x_1^{(k)})^2 &\Delta t\sum x_1^{(k)}x_2^{(k)}\\ -\Delta t\sum x_2^{(k)} &\Delta t\sum x_1^{(k)}x_2^{(k)} & \Delta t\sum (x_2^{(k)})^2 \\ \end{vmatrix}, \end{eqnarray*} $

其中$\sum\text{代表}\sum\limits_{k=1}^{n}$, 矩阵$B$各个元素的代数余子式分别为$B_{11}, B_{12}, B_{13}, B_{21}, B_{22}, B_{23}, B_{31}, B_{32}, B_{33}, $

$ \begin{eqnarray*} &(A^{T}A)^{-1}= \begin{pmatrix} n\Delta t &-\Delta t\sum x_1^{(k)} &-\Delta t\sum x_2^{(k)}\\ -\Delta t\sum x_1^{(k)} &\Delta t\sum (x_1^{(k)})^2 &\Delta t\sum x_1^{(k)}x_2^{(k)}\\ -\Delta t\sum x_2^{(k)} &\Delta t\sum x_1^{(k)}x_2^{(k)} & \Delta t\sum (x_2^{(k)})^2 \\ \end{pmatrix}^{-1} =\begin{pmatrix} \frac{B_{11}}{B} &\frac{B_{21}}{B} &\frac{B_{31}}{B} \\ \frac{B_{12}}{B} &\frac{B_{22}}{B} &\frac{B_{32}}{B} \\ \frac{B_{13}}{B} &\frac{B_{23}}{B} &\frac{B_{33}}{B}\\ \end{pmatrix}, \\ & A^{T}y_1 = \begin{pmatrix} \sum y_1^{(k)}\sqrt{\Delta t}\\ -\sum x_1^{(k)}y_1^{(k)}\sqrt{\Delta t}\\ -\sum x_2^{(k)}y_1^{(k)}\sqrt{\Delta t}\\ \end{pmatrix}, \quad A^{T}y_2 = \begin{pmatrix} \sum y_2^{(k)}\sqrt{\Delta t}\\ -\sum x_1^{(k)}y_2^{(k)}\sqrt{\Delta t}\\ -\sum x_2^{(k)}y_2^{(k)}\sqrt{\Delta t}\\ \end{pmatrix}. \end{eqnarray*} $

故(3.1)式可进一步化简为

$ \begin{eqnarray*} & (A^{T}A)^{-1}(A^{T}y_1) = \begin{pmatrix} \frac{B_{11}}{B}\sum y_1^{(k)}\sqrt{\Delta t}-\frac{B_{21}}{B}\sum x_1^{(k)}y_1^{(k)}\sqrt{\Delta t}-\frac{B_{31}}{B}\sum x_2^{(k)}y_1^{(k)}\sqrt{\Delta t} \\ \frac{B_{12}}{B}\sum y_1^{(k)}\sqrt{\Delta t}-\frac{B_{22}}{B}\sum x_1^{(k)}y_1^{(k)}\sqrt{\Delta t}-\frac{B_{32}}{B}\sum x_2^{(k)}y_1^{(k)}\sqrt{\Delta t} \\ \frac{B_{13}}{B}\sum y_1^{(k)}\sqrt{\Delta t}-\frac{B_{23}}{B}\sum x_1^{(k)}y_1^{(k)}\sqrt{\Delta t}-\frac{B_{33}}{B}\sum x_2^{(k)}y_1^{(k)}\sqrt{\Delta t} \\ \end{pmatrix}, \\ & (A^{T}A)^{-1}(A^{T}y_2) = \begin{pmatrix} \frac{B_{11}}{B}\sum y_2^{(k)}\sqrt{\Delta t}-\frac{B_{21}}{B}\sum x_1^{(k)}y_2^{(k)}\sqrt{\Delta t}-\frac{B_{31}}{B}\sum x_2^{(k)}y_2^{(k)}\sqrt{\Delta t} \\ \frac{B_{12}}{B}\sum y_2^{(k)}\sqrt{\Delta t}-\frac{B_{22}}{B}\sum x_1^{(k)}y_2^{(k)}\sqrt{\Delta t}-\frac{B_{32}}{B}\sum x_2^{(k)}y_2^{(k)}\sqrt{\Delta t} \\ \frac{B_{13}}{B}\sum y_2^{(k)}\sqrt{\Delta t}-\frac{B_{23}}{B}\sum x_1^{(k)}y_2^{(k)}\sqrt{\Delta t}-\frac{B_{33}}{B}\sum x_2^{(k)}y_2^{(k)}\sqrt{\Delta t} \\ \end{pmatrix}. \end{eqnarray*} $

为了获得参数的置信区间, 需要估计参数$\beta\text{和}\eta$的方差, 方差估计公式为

$ {\rm{Var}}(\hat{\beta})=(A^TA)^{-1}\sigma_1^2, \quad\quad {\rm{Var}}(\hat{\eta})=(A^TA)^{-1}\sigma_2^2. $ (3.2)

参数$\sigma_1^2, \sigma_2^2$由残差均方公式

$ \begin{eqnarray} \hat{\sigma}^2=\frac{(Y-X\hat{\theta})^T(Y-X\hat{\theta})}{n-p} \end{eqnarray} $ (3.3)

估计得到, $p$为被估计参数个数.由(3.1), (3.3)式可得

$ \hat{\sigma}_1^2=\dfrac{y_1^Ty_1-y_1^TA\hat{\beta}}{n-3}, \quad \text{其中}\, \hat{\beta}=(A^{T}A)^{-1}(A^{T}y_1), $ (3.4)
$ \hat{\sigma}_2^2=\dfrac{y_2^Ty_2-y_2^TA\hat{\eta}}{n-3}, \quad \text{其中}\, \hat{\eta}=(A^{T}A)^{-1}(A^{T}y_2). $ (3.5)

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

$ \left\{ \begin{aligned} \hat{\sigma}_1^2=&\frac{1}{n-3}\Big[\sum(y_1^{(k)})^2-\frac{\Delta t}{B} \big( \sum\limits_{j=1}^3B_{1j}(\sum y_1^{(k)})^2 +\sum\limits_{j=1}^3B_{2j}(\sum x_1^{(k)}y_1^{(k)})^2 \\ & +\sum\limits_{j=1}^3B_{3j}(\sum x_2^{(k)}y_1^{(k)})^2\big) \Big], \\ \hat{\sigma}_2^2=&\frac{1}{n-3}\Big[\sum(y_2^{(k)})^2-\frac{\Delta t}{B} \big( \sum\limits_{j=1}^3B_{1j}(\sum y_2^{(k)})^2 +\sum\limits_{j=1}^3B_{2j}(\sum x_1^{(k)}y_2^{(k)})^2 \\ & +\sum\limits_{j=1}^3B_{3j}(\sum x_2^{(k)}y_2^{(k)})^2\big) \Big]. \end{aligned} \right. $ (3.6)

定理3.1  方程组(3.6)中$\hat{\sigma}_1^2, \hat{\sigma}_2^2$分别为方程组(2.1)中参数$\sigma_1^2, \sigma_2^2$的渐进无偏估计, 即当$n\rightarrow \infty, $

$ \hat{\sigma}_1^2 \rightarrow\sigma_1^2 \, \, a.s \, , \quad \hat{\sigma}_2^2 \rightarrow\sigma_2^2 \, \, a.s. $

 类似文献[7]证明.

利用$\hat{\sigma}_1^2, \hat{\sigma}_2^2$分别估计(3.2)中参数$\sigma_1^2, \sigma_2^2$, 得到被估参数的方差

$ \begin{eqnarray*} {\rm{Var}}(\hat{\beta})=\begin{pmatrix} \frac{B_{11}}{B} &\frac{B_{21}}{B} &\frac{B_{31}}{B} \\ \frac{B_{12}}{B} &\frac{B_{22}}{B} &\frac{B_{32}}{B} \\ \frac{B_{13}}{B} &\frac{B_{23}}{B} &\frac{B_{33}}{B}\\ \end{pmatrix}\hat{\sigma}_1^2 \, , \quad\quad {\rm{Var}}(\hat{\eta})=\begin{pmatrix} \frac{B_{11}}{B} &\frac{B_{21}}{B} &\frac{B_{31}}{B} \\ \frac{B_{12}}{B} &\frac{B_{22}}{B} &\frac{B_{32}}{B} \\ \frac{B_{13}}{B} &\frac{B_{23}}{B} &\frac{B_{33}}{B}\\ \end{pmatrix}\hat{\sigma}_2^2. \end{eqnarray*} $
4 区间估计

如果$\sigma^2$已知, 由最小乘法回归理论[8], 参数估计值$\hat{\beta}, \hat{\eta}$的各分量都满足正态分布.当观察值的数量$n$足够大时, 可由$\hat{\sigma}^2$代替$\sigma^2, $因此$(1-\alpha)$的置信区间(CIs)分别为

$ \begin{aligned} &\hat{r}_1 \pm z_{\frac{\alpha}{2}}\sqrt{{\rm{Var}}(\hat{r}_1)}=\frac{\sqrt{\Delta t}}{B} \big(B_{11}\sum y_1^{(k)}-B_{21}\sum x_1^{(k)}y_1^{(k)}-B_{31}\sum x_2^{(k)} y_1^{(k)} \big) \pm z_{\frac{\alpha}{2}}\sqrt{\frac{B_{11}}{B}\hat{\sigma}_1^2}, \\ & \hat{a}_{11} \pm z_{\frac{\alpha}{2}}\sqrt{{\rm{Var}}(\hat{a}_{11})}=\frac{\sqrt{\Delta t}}{B} \big(B_{12}\sum y_1^{(k)}-B_{22}\sum x_1^{(k)}y_1^{(k)}-B_{32}\sum x_2^{(k)} y_1^{(k)} \big) \pm z_{\frac{\alpha}{2}}\sqrt{\frac{B_{22}}{B}\hat{\sigma}_1^2}, \\ & \hat{a}_{12} \pm z_{\frac{\alpha}{2}}\sqrt{{\rm{Var}}(\hat{a}_{12})}=\frac{\sqrt{\Delta t}}{B} \big(B_{13}\sum y_1^{(k)}-B_{23}\sum x_1^{(k)}y_1^{(k)}-B_{33}\sum x_2^{(k)} y_1^{(k)} \big) \pm z_{\frac{\alpha}{2}}\sqrt{\frac{B_{33}}{B}\hat{\sigma}_1^2}, \\ & \hat{r}_2 \pm z_{\frac{\alpha}{2}}\sqrt{{\rm{Var}}(\hat{r}_2)}=\frac{\sqrt{\Delta t}}{B} \big(B_{11}\sum y_2^{(k)}-B_{21}\sum x_1^{(k)}y_2^{(k)}-B_{31}\sum x_2^{(k)} y_2^{(k)} \big) \pm z_{\frac{\alpha}{2}}\sqrt{\frac{B_{11}}{B}\hat{\sigma}_2^2}, \\ & \hat{a}_{21} \pm z_{\frac{\alpha}{2}}\sqrt{{\rm{Var}}(\hat{a}_{21})}=\frac{\sqrt{\Delta t}}{B} \big(B_{12}\sum y_2^{(k)}-B_{22}\sum x_1^{(k)}y_2^{(k)}-B_{32}\sum x_2^{(k)} y_2^{(k)} \big) \pm z_{\frac{\alpha}{2}}\sqrt{\frac{B_{22}}{B}\hat{\sigma}_2^2}, \\ & \hat{a}_{22} \pm z_{\frac{\alpha}{2}}\sqrt{{\rm{Var}}(\hat{a}_{22})}=\frac{\sqrt{\Delta t}}{B} \big(B_{13}\sum y_2^{(k)}-B_{23}\sum x_1^{(k)}y_2^{(k)}-B_{33}\sum x_2^{(k)} y_3^{(k)} \big) \pm z_{\frac{\alpha}{2}}\sqrt{\frac{B_{33}}{B}\hat{\sigma}_2^2}, \\ \end{aligned} $

其中$z_{\frac{\alpha}{2}}$是上$\frac{\alpha}{2}$值的标准正态随机变量, 例如: $z_{0.025}=1.96$为95%的置信区间.

下面给出具体数据, 进行参数的估计.

例1  在模型(1.2)式中, 我们选择满足种群持久生存的数据进行模拟, 即两种群随机持久存在的情形, 其它状态的数据也可获得类似结论.假设参数真实值分别为$r_1=1, a_{11}=0.3, a_{12}=0.2, \sigma_1=0.04, r_2=1.2, a_{21}=0.3, a_{22}=0.4, \sigma_2=0.05, $初始值为$x_1(0)=1.5, x_2(0)=2.$应用这些参数的值, 通过Euler-Maruyama方法离散, 我们获得三组模拟数值$x_1^{k}(t), x_2^{k}(t), $数量分别为$(A)n=5000; (B)n=10000; (C)n=50000; $保存这些数据为模拟数据集.对每一组数据, 给步长$ \Delta t=0.02, 0.05, 0.1$进行模拟, 给出一些模拟结果, 以便比较真实值与估计值.在表 1-4中, 样本的大小以符号“Size n”表示, 且在表中第一列给出.其中表 1, 2样本大小从5000增加到50000,在每一组模拟数据中, 分别以三种不同步长进行模拟.模拟数据表明参数的估计值与真实值没有明显区别, 绝对误差与样本大小有关, 与步长$\Delta t$无关, 且随着样本数量的增大, 绝对误差越来越小.

表 1 $(r_1, a_{11}, a_{12})$的估计模拟结果(真实值: $r_1=1, a_{11}=0.3, a_{12}=0.2$)

表 2 $(r_2, a_{21}, a_{22})$的估计模拟结果(真实值: $r_2=1, a_{21}=0.3, a_{22}=0.2$)

表 3 置信水平为0.95的参数$r_1, a_{11}, a_{12}$置信区间模拟结果(真实值: $r_1=1, a_{11}=0.3, a_{12}=0.2$)

表 4 置信水平为0.95的参数$r_2, a_{21}, a_{22}$置信区间模拟结果(真实值: $r_2=1, a_{21}=0.3, a_{22}=0.2$)

其次给出各个参数置信水平为0.95的置信区间的一些模拟结果. 表 3, 4中, 给出参数的估计区间及区间长度, 模拟数据表明随着样本量的增大, 置信水平长度越来越小.

5 结论

本文对随机两种群Lotka-Volterra竞争模型应用最小二乘法理论进行参数估计, 获得参数$r_1, a_{11}, a_{12}, r_2, a_{21}, a_{22}, \sigma_1, \sigma_2$的估计值及估计区间.通过观察结果, 得到影响置信区间长度的主要因素.结果表明, 参数估计区间的长度随着样本数量n的增加而减小, 不依赖步长$\Delta t$的长度.例1给出了具体的数值模拟.极大似然估计与贝叶斯估计方法是另外两种参数估计的方法, 在将来的工作中, 将应用这些方法到随机两种群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. DOI:10.1016/j.automatica.2003.10.001
[3] Timmer J. Parameter estimation in nonlinear stochastic differential equatuons[J]. Chaos Sol. Fract., 2000, 11(15): 2571–2578. DOI:10.1016/S0960-0779(00)00015-1
[4] Zhu C, Yin G. On competitive Lotka-Volterra model in random environments[J]. J. Math. Anal. Appl., 2009, 357(1): 154–170. DOI:10.1016/j.jmaa.2009.03.066
[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. DOI:10.1016/j.jmaa.2011.12.049
[6] 王佳, 丁洁丽. Logistic回归模型中参数极大似然估计的二次下界算法及其应用[J]. 数学杂志, 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. DOI:10.1007/s11203-014-9091-8
[8] Rawlings JO. Applied regression analysis:a research tool[M]. Belmont, CA: Wadsworth, 1988.
[9] 洪志敏, 闫在在. Volterra积分方程的蒙特卡罗数值求解方法[J]. 数学杂志, 2016, 36(2): 425–436.