Microsoft Word - Stat040225-arangEq-1.doc SQU Journal For Science, 9 (2004) 53-66 © 2004 Sultan Qaboos University 53 Recursive Estimation of a Discrete- Time Lotka-Volterra Predator-Prey Model Lakhdar Aggoun Department of Mathematics and Statistics, Sultan Qaboos University, P. O. Box 36, Alkhod 123, Sultanate of Oman, Email: laggoun@squ.edu.om. فلتره للصائدة والمصطادة-التقدير الشرطي المتكرر لنموذج لوتكه األخضر عقون هذا المقال يعالج مسألة مهمة في اإلحصاء البيني وهي مسألة حساب التعداد الحيواني لنـوعين مـن : خالصة طلق سراحها من كل نوعية تؤخذ عينة عشوائية وتميز بإشارة خاصة ثم ي . الحيوانات وهي الصائدة و المصطادة وبعد فترة تؤخذ عينة عشوائية من القطيعين ويحصي منها األفراد المميزين باإلشارة الخاصة . داخل بقية القطيع ومن هذه المعلومات تستعمل طرق القياس لتقدير التوزيع الشرطي المتكـرر لمجموعـة الحيوانـات الصـائدة .والمصطادة ABSTRACT: In this paper, using hidden Markov models, we estimate the number of individuals in a two-species (predator-prey) animal population using partial information provided by the so-called capture-recapture technique. Random samples of individuals are captured, tagged in some way and released. After some time other random samples are taken and the marked individuals are observed. Using this information, we estimate (recursively) the sizes of the two populations. Also, using the Expectation Maximization (EM) algorithm, the parameters of the model are updated. KEYWORDS: Hidden markov models, predator-prey, capture-recapture, EM algorithm. AMS subject classification: 60J27, 93E11. 1. Introduction ne of the first models to incorporate interactions between predators and preys was proposed in 1925 by the American biophysicist Alfred Lotka and the Italian mathematician Vito Volterra . Vito Volterra (1860-1940) was a famous Italian mathematician who retired from a distinguished career in pure mathematics in the early 1920s. His son-in-law, Humberto D’Ancona, was a biologist who studied the populations of various species of fish in the Adriatic. Volterra developed a series of models for interactions of two or more species. O LAKHDAR AGGOUN 54 Alfred J. Lotka (1880-1949) was an American mathematical biologist (and later actuary) who formulated many of the same models as Volterra, independently and at about the same time. His primary example of a predator-prey system comprised a plant population and an herbivorous animal dependent on that plant for food. Here we consider a discrete-time stochastic version of the Lotka-Volterra model. Hidden Markov models (Elliott et al. 1995) have been used extensively in many areas of science and technology. In this paper we are extending the use of this powerful tools (see Aggoun and Elliott 1998 for a single species model) to estimate the hidden number of individuals in a multi-species animal population using partial information provided by the so-called capture-recapture technique (see Seber 1982, for instance). Two random samples of individuals one from the predator population and one from the prey population are captured, tagged or marked in some way, and then released. After allowing time for the marked and unmarked to mix sufficiently, two second simple random samples from both populations are taken and the marked ones are observed. At epoch write Z for the prey population size, z for the number of marked and released preys, 1 k k zz == ∑ for the total number of captured and marked preys up to time k , R for the sample size, z for the number of available marked individuals for sampling and Zy for the number of captured (or recaptured) marked individuals. Similarly write X for the predator population size, x for the number of marked and released predators, 1 k k xx == ∑ for the total number of captured and marked predators up to time k , F for the sample size, x for the number of available marked predators for sampling and Xy for the number of captured (or recaptured) marked predators. 2. Model Assumptions and Recursive Estimation All random variables are defined initially on a probability space ( )F PΩ, , . All the filtrations defined here are assumed to be complete. Write ( )Z XkG Z z y R X x y F kσ= , , , , , , , ≤ , and ( , ) Z X kY y y kσ= , ≤ . kX and kZ represent the number of predators and preys, respectively, that are alive at time period k, then a (discrete time) Lotka-Volterra type model is: 1 1 1k k k k k kX X aZ X bX vσ+ += + − + (1) 1 2 1k k k k k kZ Z cZ dZ X wσ+ += + − + , where the parameters are defined by: a is the efficiency of turning predated preys into predators. b is the natural death rate of predators in the absence of food (preys), c is the natural growth rate of preys in the absence of predation, d is the death rate per encounter of preys due to predation, v and w are sequences of independent random variables with some (either discrete or continuous with finite supports) densities kφ and kϑ respectively and 1σ , 2σ are some positive real numbers. The random RECURSIVE ESTIMATION OF A DISCRETE-TIME 55 variables v and w indicate other sources of variations in the populations like death caused by old age, diseases, births etc. It is assumed here that kx and kz are random variables with binomial distributions with parameters ( )Xk px , and ( )Zk pz , respectively. The observed random variables Z Xk ky y, are assumed to have joint distribution: 1( ) X Z k k k k k k k k kP y m y n G z Z R x X F−= , = | , , , , , , ( ) (1 ) ( ) (1 )k kk k F m R nm nk k k k k k k k F R x x z z m n X X Z Z − −  = − −      (5) Let 0 1λ = . For 1≥ and for suitable density functions χ , v write X Zλ λ λ= , where 1 ( ) 1 (1 ) ( ) (1 ) ( ) 2 X Xxx y y FX x X XF x X x x p p v X X σ χ λ φ − +− − − + = − − (6) 2 ( ) 1 (1 ) ( ) (1 ) ( ) 2 Z Zzz y y RZ z Z ZR z v Z z z p p w Z Z σ λ ϑ − +− − − + = − − (7) and 0 k k λ=Λ =∏ . Lemma 1. The process kΛ is a G -martingale. Proof. 1 1 1[ ] [ ] X Z k k k k k kE G E Gλ λ− − −Λ | = Λ | , so we must show that 1[ ] 1 X Z k k kE Gλ λ −| = . However using repeated conditioning it is enough to show that 1[ ] 1 Z k kE Gλ −| = , say. 2 1 1 2 1 1 2 ( ) 1 [ ] [ ( ) (1 ) ( ) 2 (1 ) ] ( )1 [ (1 ) [( ) ( )2 (1 ) ] ] ( )1 [ ( )2 Z Z k k k k k kk k Z kk k k k k Z k k k k y y RZ k k k k k k R z k k k k zz z Z kZ zz yk k kz Z ZR z k k k y Rk k k k k k k zk k ZR z k k v Z z z E G E w Z Z p p G v Z z E p p E w Z z G Z z R G Z v Z E p w σ λ ϑ σ ϑ σ ϑ − − − + − +− − − +− − + − − − − + | = − × − | = − × − | , , , | = 1 0 (1 ) ] k kk k R z kz Z k m R p G m − + − =   − |    ∑ LAKHDAR AGGOUN 56 2 1 2 1 0 2 1 2 1 2 1 2 ( )1 [ (1 ) ] ( )2 ( )1 [ ] ( )2 ( ) [ ] ( ) ( ) ( ) ( ) 1 ( ) kk k k k k zzk k z Z Z kz k k z k k k kz ik k k k k k k k k k k k k v Z E p p G w v Z zE G w i v Z w E G w v Z w w dv v u du w σ ϑ σ ϑ σ σ ϑ σ σ ϑ ϑ − +− − − = − − − = − |   = |    + = | + = = = ∑ ∫ ∫ A new probability measure Q can be defined by setting | k k dQ E G dP   = Λ   . The point here is that: Lemma 2. Under the new probability measure Q , kX , kx , X ky kZ , kz and Z ky are sequences of independent random variables which are independent of each other. Further, kX has density kχ , kx has distribution 1 ( ) 2k bin x , , x ky has distribution 1 ( ) 2k bin F , , kZ has density kv , kz has distribution 1 ( ) 2k bin z , and Z ky has distribution 1 ( ) 2k bin R , . Proof. We shall check the claim for the three processes kZ , kz and Z ky . For any integrable real-valued functions f , g and h and using a version of Bayes’ theorem (see Elliott et al. 1995) we can write: 1 1 1 1 2 1 [ ( ) ( ) ( ) ] [ ( ) ( ) ( ) ] [ ] [ ( ) ( ) ( ) ] ( ) 1 [ ( ) ( ) ( ) (1 ) ( ) 2 1[ ( ) ( ) (1 ) ] 2 kk k k Z Z k k k k Z Z k k k k k Q k k k k k k Z k k k k k zzZ k k n k k k Z Zn k k y y Rk kZ k k k kk R k k E f Z g z h y G E f Z g z h y G E G E f Z g z h y G v Z E f Z g z h y p p w z zE h y G Z z R G Z Z λ σ ϑ − − − − − +− − − − Λ | | = Λ | = | = − − | , , , | 1 2 1 0 ] ( ) 1 [ ( ) ( ) ( ) (1 ) ( ) 2 1[ ( )( ) (1 ) ( ) (1 ) ] ] 2 kk k k k k k k k zzZ k k n k k k Z Zn k k R km R m Rk k k km m kR m k k k k v Z E f Z g z h y p p w Rz z z zh m G mZ Z Z Z σ ϑ − − +− − + −− − = = −   − − |     ∑ RECURSIVE ESTIMATION OF A DISCRETE-TIME 57 2 1 0 2 1 1 1 1 ( ) 1 [ ( ) ( ) ( ) (1 ) ( ) 2 1[ ( ) ] ] 2 ( ) 1 [ ( )] [ ( ) ( ) (1 ) ] ( ) 2 [ ( )] [ ( )] [ ( kk k k k k kk k k zzZ k k n k k k Z Zn k k R k kR m zzZ k k n Q k k k Z Z kn k k Z Q k Q k k k k v Z E f Z g z h y p p w R h m G m v Z E h y E f Z g z p p G w E h y E g z E f Z cZ dZ σ ϑ σ ϑ − +− − = − +− − − − − = −   × |     = − | = + − ∑ 2 2 1 1 1 2 1 1 1 1 2 2 1 1 1 2 ) ( ) ] ( ) [ ( )] [ ( )] ( ) ( ) [ ( )] [ ( )] ( ) ( ) [ ( )] [ ( )] [ ( )], k k k k k k k k k k k Z Q k Q k k k k k k k k k k Z Q k Q k k Z Q k Q k Q k X w v Z cZ dZ X w G w E h y E g z f Z cZ dZ X w v Z cZ dZ X w dw E h y E g z f u v u du E h y E g z E f Z σ σ σ ϑ σ σ σ − − − − − − − − − − + + − + × | = × + − + × + − + = = ∫ ∫ where k R m        = ( ) ! ! ! k k R m R m− . That is, under Q the three processes are independent sequences of random variables with the desired distributions. Using this fact we derive a recursive equation for the unnormalized conditional distribution of kZ and kX given kY . For any measurable test function f consider: 1 1 [ ( ) ] [ ( ) ] [ ] Q k k k k k k k Q k k E f X Z Y E f X Z Y E Y − − , Λ | , | = . Λ | (8) The denominator of (8) being a normalizing factor we focus only on the expectation under Q in the numerator. Write 1[ ( ) ] ( ) ( )Q k k k k kE f X Z Y f x z q x z dxdz −, Λ | = , , .∫ (9) Theorem 1. The unnormalized conditional joint probability density function of the populations’ sizes given by the dynamics in (1) and (2) follows the recursion: LAKHDAR AGGOUN 58 1( ) ( ) ( ) ( )k k k kq x z A x z B u v x z q u v dudv−, = , , , , ,∫ (10) where 0 01 2 2 ( ) (1 ) (1 ) k k k k k k R F x z i ji jz xk k k Z Z X X j i z xA x z p p p p i jσ σ + − − = =    , = − −      ∑∑ ( ) (1 ) ( ) (1 ) X X Z Z k k k k k ky F y y R yj j i i x x z z − −× − − (12) 2 1 ( ) ( ) ( )k k k z v cv dvu x u avu bu B u v x z ϑ φ σ σ − − + − − + , , , = (13) (Note we take 00 1= .) Proof. In view of Lemma 2 the left hand side of (9) is: 1 1 1 0 0 1 1 1 1 2 2 [ ( ) ] 1 1 2 2 2 [ (1 ) (1 ) ( )( ) (1 ) ( ) (1 ) ( ) ( ) ( ) ( k kk k k k k k k k X X Z Z k k k k k k Q k k k k k R Fz x x z x z i ji jz xk k Q Z Z X X j i y F y y R y k k k k k k k k E f X Z Y z xE p p p p i j j j i if x z x x z z z Z cZ dZ X v z v z x X λ ϑ σ σ φ − − − + + + − − = = − − − − − − = , Λ | =    × − −      × , − − − − + − × ∑∑ ∫ 1 1 1 1 1 1 1 1 0 01 2 ) ( ) ( ) (1 ) (1 ) ] 2 (1 ) (1 ) ( )( ) (1 ) ( ) (1 ) ( k k k k k k k k X X Z Z k k k k k k k k k k k k i ji jz x Z Z X X k k R F x z i ji jz xk k Z Z X X j i y F y y R y k aZ X bX x x x p p p p dxdz Y z x p p p p i j j j i if x z x x z z z v cv d σ σ χ σ σ ϑ − − − − − − − − + − − = = − − − + × − − Λ |    = − −      × , − − − − + ∑∑ ∫ ∫ 1 2 1 ) ( ) ( )k k vu x u avu bu q u v dudvdxdzφ σ σ − − − + , Comparing this last expression with the right hand side of (9) gives the result. RECURSIVE ESTIMATION OF A DISCRETE-TIME 59 Remark 1. The normalized conditional joint density function of ( )k kX Z, is simply ( ) ( ) k k q x z q u v dudv , ,∫ . The initial (normalized) probability density of 0 0( )X Z, , prior to sampling, is 0 ( )π . , so 0 0( ) ( )q x z x zπ, = , . 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 01 2 1 1 0 2 1 2 ( ) (1 ) (1 ) ( ) (1 ) ( ) (1 ) ( ) ( ) ( ) . X X Z Z R F x z jii j xz Z Z X X j i y F y y R y xzq x z p p p p i j j j i i x x z z z v cv dvu x u avu bu u v dudv σ σ ϑ φ π σ σ + −− = = − −    , = − −      × − − − − + − − + . ∑∑ ∫ And further estimates follow from (10). If the distribution of 0 0( )X Z, is a delta function concentrated at ( )A B, say 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 01 2 1 1 2 1 2 ( ) (1 ) (1 ) ( ) (1 ) ( ) (1 ) ( ) ( ). X X Z Z R F x z jii j xz Z Z X X j i y F y y R y xzq x z p p p p i j j j i i x x z z z B cB dBA x A aBA bA σ σ ϑ φ σ σ + −− = = − −    , = − −      × − − − − + − − + ∑∑ 3. Parameter Revision Here we shall assume, for simplicity, that the random variables v and w in our model are standard normal (means 0 and standard deviations 1). The EM algorithm, (Baum and Petrie 1966, Dempster et al. 1977) is a widely used iterative numerical method for computing maximum likelihood parameter estimates of partially observed models such as linear Gaussian state space models. For such models, direct computation of the MLE is difficult. The EM algorithm has the appealing property that successive iterations yield parameter estimates with non- decreasing values of the likelihood function. Suppose that we have observations 1 Ky … y, , available, where K is a fixed positive integer. Let { }Pθ θ, ∈Θ be a family of probability measures on ( )FΩ, , all absolutely continuous with respect to a fixed probability measure 0P . The log-likelihood function for computing an estimate of the parameter θ based on the information available in KY is 0 0 ( ) logK K dP L E Y dP θθ            = | , and the maximum likelihood estimate (MLE) is defined by ˆ argmax ( )KLθθ θ∈Θ∈ . LAKHDAR AGGOUN 60 Let 0θ̂ be the initial parameter estimate. The EM algorithm generates a sequence of parameter estimates { }jθ , 1j ≥ , as follows. Each iteration of the algorithm consists of two steps: Step 1. (E-step). Set ˆ jθ θ= and compute ( )Q θ θ, , where ( ) log K dP Q E Y dP θ θ θ θ θ            , = | . Step 2. (M-step). Find 1ˆ argmax ( )jj Qθ θ θθ ∈Θ+ ∈ , . Using Jensen’s inequality it can be shown (see Theorem 1 in Dempster et al. 1977) that the sequence of model estimates ˆ{ jθ , 1}j ≥ from the EM algorithm are such that the sequence of likelihoods ˆ{ ( )}K jL θ , 1j ≥ is monotonically increasing with equality if and only if 1ˆ ˆj jθ θ+ = . Sufficient conditions for convergence of the EM algorithm are given in Wu (1983). We briefly summarize them here: 1. The parameter space Θ is a subset of some finite dimensional Euclidean space rR . 2. 0ˆ{ ( ) ( )}K KL Lθ θ θ∈Θ: ≥ is compact for any 0ˆ( )KL θ > −∞ 3. KL is continuous in Θ and differentiable in the interior of Θ . (As a consequence of (1), (2) and (iii), clearly ˆ( )K jL θ is bounded from above). 4. The function ˆ( )jQ θ θ, is continuous in both θ and ˆ jθ . Then by Theorem 2 in Wu (1983), the limit of the sequence EM estimates ˆ{ }jθ has a stationary point θ of KL . Also ˆ{ ( )}K jL θ converges monotonically to ( )t tLL θ= for some stationary point θ . To make sure that tL is a maximum value of the likelihood, it is necessary to try different initial values 0θ̂ . The model in (1) and (2) is determined by the parameters a , b , c and d which need to be updated as new information is obtained. These parameters are estimated using the expectation maximization (EM) algorithm. Maximum likelihood estimation of the parameters via the EM algorithm requires computation of the filtered estimates of quantities such as ( ) 2 2 1 11 kI kT X Z− −== ∑ , ( ) 2 1 11 kII kT X Z− −== ∑ , ( ) 1 11 kIII kT X X Z− −== ∑ , ( ) 2 11 kIV kT X −== ∑ , ( ) 11 kV kT X X −== ∑ , ( ) 2 1 11 kI kS X Z− −== ∑ , RECURSIVE ESTIMATION OF A DISCRETE-TIME 61 ( ) 1 11 kII kS X Z Z− −== ,∑ ( ) 2 1 kIII kS Z== ,∑ ( ) 11 kIV kS Z Z −== ∑ , 1 k k kU x== ∑ and 1 k k kV z== .∑ For M I … V= , , write 1 ( ) ( ) 1 [ ( ) ] [ ( ) ] [ ] M Q k k k k kM k k k k Q k k E T I X dx Z dz Y E T I X dx Z dz Y E Y − − Λ ∈ , ∈ | ∈ , ∈ | = Λ | , ( ) 1 ( )( ) [ ( ) ]M Mk Q k k k k kx z E T I X dx Z dz Yβ −, = Λ ∈ , ∈ | , 1 ( ) ( ) ( )( ) 1 [ ] ( ) [ ] [ ] ( ) M M MQ k k kM k k k k Q k k k E T Y x z dxdz E T Y TE Y q x z dxdz β− − Λ | , | = = Λ | , ∫ ∫ . For N I … IV= , , . 1 ( ) ( ) 1 [ ( ) ] [ ( ) ] [ ] N Q k k k k kN k k k k Q k k E S I X dx Z dz Y E S I X dx Z dz Y E Y − − Λ ∈ , ∈ | ∈ , ∈ | = Λ | , ( ) 1 ( )( ) [ ( ) ]N Nk Q k k k k kx z E S I X dx Z dz Yη −, = Λ ∈ , ∈ | , 1 ( ) ( ) ( )( ) 1 [ ] ( ) [ ] [ ] ( ) N N NQ k k kN k k k k Q k k k E S Y x z dxdz E S Y SE Y q x z dxdz η− − Λ | , | = = Λ | , ∫ ∫ . First we compute ML estimates of the parameters 1 2( )a b c dθ σ σ= , , , , , given the history in kY . Now the expression for ( )Q θ θ, is derived. To update the set of parameters from θ to θ , the following density is introduced 1 k k G dP dP θ θ γ = | = ,∏ 1 1 1 1 1 1 1 1 1 2 1 2 1 1 1 1 1 1 1 1 1 2 1 2 ( ) ( ) . ( ) ( ) k k k k k k k k k k k k k k k k k k k k X X aZ X bX Z Z cZ dZ X X X aZ X bX Z Z cZ dZ X σ φ σ ϑ σ σ γ σ φ σ ϑ σ σ − − − − − − − − − − − − − − − − − − + − − + = − − + − − + Now 21 1 1 1 1 2 1 1 21 1 1 1 1 2 1 log log log ( ) 2 1 ( ) ( ) ( ) 2 k k k k k k k G k k k k k k k k k dP X X aZ X bX E Y k k E Y dP Z Z dZ X cZ E Y R Q θ θ θ θ θ σ σ σ θ θ θ σ        − − − −        =         − − − −     =   − − + | | = − − − | − + − − | + = , , ∑ ∑ LAKHDAR AGGOUN 62 where ( )R θ does not involve θ . To implement the M-step set the derivatives 0 Q θ ∂ = ∂ . This yields ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )2 2 1 1 1 1 1 ( 1) 1 ( ) 2 ( ) 1 ( III II k k I k III II V I k k k k IV I II k k k IV III I k k k III k III I IV I II k k k k k I I III k k k k k k k k k bT Ta T T T T Tb T T T dS S Sc S S S S S Sd S ST E X X aZ X bX k θ σ − − − − = + − = , − − = , − − + = , − + = − = − − +∑ 21 2 2 2 1 1 1 1 1 ) 1 ( ) k k k k k k k k Y E Z Z dZ X cZ Y k θ σ                − − − −   =  | , = − + − | .∑ For any “test” function g , write 1 ( ) ( )[ ( ) ] ( ) ( )M MQ k k k k k kE T g X Z Y x z g x z dxdz M I … Vβ −Λ , | = , , , = , , .∫ 1 ( ) ( )[ ( ) ] ( ) ( )N NQ k k k k k kE S g X Z Y x z g x z dxdz N I … IVη −Λ , | = , , , = , , .∫ Theorem 1. For 1k ≥ , the unnormalized densities ( ) ( )Mk x zβ , , ( ) ( )Nk x zη , , M I … V= , , , N I … IV= , , are given by the following recursions. ( ) ( ) 1 2 2 1 ( ) ( ) 1 2 1 ( ) ( ) 1 ( ) ( )[ ( ) ( ) ( ) ( ) ] ( ) ( )[ ( ) ( ) ( ) ( ) ] ( ) ( )[ ( ) ( I I k k k k k k II I k k k k k k III I k k k k x z A x z B u v x z u v dudv B u v x z u v q u v dudv x z A x z B u v x z u v dudv B u v x z u vq u v dudv x z A x z B u v x z u β β β β β β − − − − − , = , , , , , + , , , , , = , , , , , + , , , , , , = , , , , , ∫ ∫ ∫ ∫ ∫ )v dudv RECURSIVE ESTIMATION OF A DISCRETE-TIME 63 1 ( ) ( ) 1 2 1 ( ) ( ) 1 1 ( ) ( ) ( ) ] ( ) ( )[ ( ) ( ) ( ) ( ) ] ( ) ( )[ ( ) ( ) ( ) ( ) ] ( ) ( k k IV I k k k k k k V I k k k k k k I k k x B u v x z uvq u v dudv x z A x z B u v x z u v dudv B u v x z u q u v dudv x z A x z B u v x z u v dudv x B u v x z uq u v dudv x z A x β β β β η − − − − − + , , , , , , = , , , , , + , , , , , = , , , , , + , , , , , , = ∫ ∫ ∫ ∫ ∫ ( ) 1 2 1 ( ) ( ) 1 1 ( ) ( ) 1 2 )[ ( ) ( ) ( ) ( ) ] ( ) ( )[ ( ) ( ) ( ) ( ) ] ( ) ( )[ ( ) ( ) ( I k k k k II II k k k k k k III III k k k k k z B u v x z u v dudv B u v x z uv q u v dudv x z A x z B u v x z u v dudv z B u v x z uvq u v dudv x z A x z B u v x z u v dudv z q x η η η η η − − − − − , , , , , + , , , , , = , , , , , + , , , , , = , , , , , + , ∫ ∫ ∫ ∫ ∫ ( ) ( ) 1 1 )] ( ) ( )[ ( ) ( ) ( ) ( ) ] IV IV k k k k k k z x z A x z B u v x z u v dudv z B u v x z vq u v dudv η η − − , = , , , , , + , , , , ∫ ∫ where ( )kA x z, and ( )kB u v x z, , , are given in (12) and (13) and kq is given recursively in (10). Proof. First note that ( ) ( ) 2 21 1 1 I I k k k kT T X Z− − −= + . Therefore 1 ( ) 1 ( ) 1 2 2 1 1 1 [ ( ) ] [ ( ) ] [ ( ) ] I Q k k k k k I Q k k k k k Q k k k k k k E T g X Z Y E T g X Z Y E X Z g X Z Y − − − − − − = Λ , | = Λ , | + Λ , | . The first expectation is simply 1 ( ) 1 0 01 2 [ ( ) ] 2 (1 ) (1 ) ( )( ) (1 ) ( ) (1 ) k k k k k k X X Z Z k k k k k k I Q k k k k k R F x z i ji jz xk k Z Z X X j i y F y y R y E T g X Z Y z x p p p p i j j j i ig x z x x z z σ σ − − + − − = = − − Λ , |    = − −      × , − − ∑∑ ∫ LAKHDAR AGGOUN 64 2 1 ( ) 1 ( ) ( ) ( ) . k k I k z v cv dvu x u avu bu u v dudvdxdz ϑ φ σ σ β − − − + − − + × , ∫ And 1 2 2 1 1 0 01 2 2 1 2 2 1 [ ( ) ] 2 (1 ) (1 ) ( )( ) (1 ) ( ) (1 ) ( ) ( ) ( ) k k k k k k X X Z Z k k k k k k Q k k k k k k R F x z jii j xzk k Z Z X X j i y F y y R y k k k E X Z g X Z Y xz p p p p i j j j i ig x z x x z z z v cv dvu x u avu bu u v q u v dudvd σ σ ϑ φ σ σ − − − + −− = = − − − Λ , |    = − −      × , − − − − + − − + × , ∑∑ ∫ ∫ xdz Using (12) and (13) yields the result. Write 1 1 [ ( ) ] [ ( ) ] [ ] Q k k k k k k k k k Q k k E U I X dx Z dz Y E U I X dx Z dz Y E Y − − Λ ∈ , ∈ | ∈ , ∈ | = Λ | , 1( ) [ ( ) ]k Q k k k k kx z E U I X dx Z dz Yζ −, = Λ ∈ , ∈ | , 1 1 [ ] ( ) [ ] [ ] ( ) Q k k k k k k k Q k k k E U Y x z dxdz E U Y UE Y q x z dxdz ζ− − Λ | , | = = Λ | , ∫ ∫ , 1 1 [ ( ) ] [ ( ) ] [ ] Q k k k k k k k k k Q k k E V I X dx Z dz Y E V I X dx Z dz Y E Y − − Λ ∈ , ∈ | ∈ , ∈ | = Λ | , 1( ) [ ( ) ]k Q k k k k kx z E V I X dx Z dz Yξ −, = Λ ∈ , ∈ | , 1 1 [ ] ( ) [ ] [ ] ( ) Q k k k k k k k Q k k k E V Y x z dxdz E V Y VE Y q x z dxdz ξ− − Λ | , | = = Λ | , ∫ ∫ . To update the parameter from Xp to Xp , the following density is introduced 1 1 . 1 k kk X k X x xxk p X X G p X X dP p p dP p p − =    − | =     −    ∏ Now RECURSIVE ESTIMATION OF A DISCRETE-TIME 65 1 log log log(1 ) log(1 ) ( ) ( ) X X k X k p p G k X X Xk k p X X X dP E Y p p pxU UdP R p Q p p           =  | | = + − + − + = , , ∑ where ( )XR p does not involve Xp . To implement the M-step set the derivatives 0 X Q p ∂ = ∂ . This yields 1 k X k Up x = = . ∑ A similar argument yields: 1 k Z k Vp z = = . ∑ For any “test” function g , write 1[ ( ) ] ( ) ( )Q k k k k k kE U g X Z Y x z g x z dxdzζ −Λ , | = , , .∫ 1[ ( ) ] ( ) ( )Q k k k k k kE V g X Z Y x z g x z dxdzξ −Λ , | = , , .∫ Theorem 2. For 1k ≥ , the unnormalized densities ( )k x zζ , and ( )k x zξ , are given by the following recursions. 1 1 ( ) ( ) ( ) ( ) ( ) 2 ( ) ( ) ( ) ( ) ( ) 2 k k k k k k k k k k k k x z A x z B u v x z u v dudv x q x z x z A x z B u v x z u v dudv z q x z ζ ζ ξ ξ − − , = , , , , , + , , = , , , , , + , ∫ ∫ Proof. First note that 1k k kU U x−= + . Therefore 1 1 1 1 [ ( ) ] [ ( ) ] [ ( ) ] Q k k k k k Q k k k k k Q k k k k k E U g X Z Y E U g X Z Y E x g X Z Y − − − − = Λ , | = Λ , | + Λ , | . The first expectation is simply LAKHDAR AGGOUN 66 1 1 0 01 2 2 1 1 [ ( ) ] 2 (1 ) (1 ) ( )( ) (1 ) ( ) (1 ) ( ) ( ) ( ) k k k k k k X X Z Z k k k k k k Q k k k k k R F x z jii j xzk k Z Z X X j i y F y y R y k k k E U g X Z Y xz p p p p i j j j i ig x z x x z z z v cv dvu x u avu bu u v dudvdxdz σ σ ϑ φ σ σ ζ − − + −− = = − − − Λ , |    = − −      × , − − − − + − − + × , . ∑∑ ∫ ∫ from which the result follows. Remark. The above theorems do not require v and w to be Gaussian. 4. References AGGOUN, L. and ELLIOTT, R.J. 1998. Recursive Estimation in Capture-Recapture Methods, Journal of Science and Technology, Sultan Qaboos University 3: 67–75. BAUM, L.E. and PETRIE, T. 1966. Statistical Inference for Probabilistic Functions of Finite State Markov Chains, Annals of the Institute of Statistical Mathematics 37: 1554-1563. DEMPSTER, A.P., LAIRD, N.M. and RUBIN, D.B. 1977. Maximum Likelihood from Incomplete Data via the EM Algorithm, Jour. of the Royal Statistical Society, B, 39: 1-38. ELLIOTT, R.J., AGGOUN, L. and MOORE, J.B. 1995. Hidden Markov Models : Estimation and Control. Applications of Mathematics. Vol. 29 Springer-Verlag, Berlin-Heidelberg-New York. LOTKA, A., J. 1925. Elements of physical biology. Baltimore: Williams & Wilkins Co. SEBER, G.A.F. 1982. The Estimation of Animal Abundance and Related Parameters. Second Edition. Edward Arnold. VOLTERRA, V. 1926. Variazioni e fluttuazioni del numero d’individui in specie animali conviventi. Mem. R. Accad. Naz. dei Lincei. Ser. VI, vol. 2. WU, C.F.J. 1983. On the Convergence Properties of the EM algorithm, The Annals of Statistics 11, pp. 95- 103. Received 25 February 2004 Accepted 26 September 2004