International audienceWe consider a Monte Carlo Markov chain (MCMC) algorithm for fisheries stock assess- ment. The biomass of this stock at a given year could be modeled as a nonlinear function of the biomass and catch for the two previous years, of different parameters (recruitment, growth rate, nat- ural mortality rate). Given a time series of annual catch and effort data, we would like to achieve the best fitting between the data and a class of non linear/non Gaussian state-space models

Campillo, Fabien

Rakotozafy, Rivo

INRIA a CCSD electronic archive server

HAL Id: hal-00652092https://hal.inria.fr/hal-00652092Submitted on 14 Dec 2011HAL is a multi-disciplinary open accessarchive for the deposit and dissemination of sci-entific research documents, whether they are pub-lished or not. The documents may come fromteaching and research institutions in France orabroad, or from public or private research centers.L’archive ouverte pluridisciplinaire HAL, estdestinée au dépôt et à la diffusion de documentsscientifiques de niveau recherche, publiés ou non,émanant des établissements d’enseignement et derecherche français ou étrangers, des laboratoirespublics ou privés.MCMC for non linear/non Gaussian state-space models:Application to fishery stock assessmentFabien Campillo, Rivo RakotozafyTo cite this version:Fabien Campillo, Rivo Rakotozafy. MCMC for non linear/non Gaussian state-space models: Applica-tion to fishery stock assessment. 7th African Conference on Research in Computer Science, Nov 2004,Hammamet, Tunisia. ￿hal-00652092￿RubriqueMCMC for non linear/non gaussianstate–space modelsApplication to fishery stock assessmentFabien Campillo* — Rivo Rakotozafy*** INRIA/IRISACampus de Beaulieu35042 RennesFRANCEFabien.Campillo@inria.fr** Faculté des SciencesUniversité de Fianarantsoa301 FianarantsoaMADAGASCARrrakotoz@syfed.refer.mgABSTRACT.We consider a Monte Carlo Markov chain (MCMC) algorithm for fisheries stock assess-ment. The biomass of this stock at a given year could be modeled as a nonlinear function of thebiomass and catch for the two previous years, of different parameters (recruitment, growth rate, nat-ural mortality rate). Given a time series of annual catch and effort data, we would like to achieve thebest fitting between the data and a class of non linear/non Gaussian state–space models.RÉSUMÉ. On applique un algorithme de Monte Carlo par Chaîne de Markov (MCMC) à l’évaluationd’un stock de pêcherie. La biomasse de ce stock à une année donnée peut être modélisée commeune fonction non linéaire des biomasses et des prises aux deux années précédentes, ainsi que dedifférents paramètres (recrutement, taux de croissance, taux de mortalité naturelle). Etant donnée unesérie temporelle de captures et d’efforts, on souhaite parvenir au meilleur accord entre ces donnéeset une classe de modèles à espace d’état non linéaire/non gaussien.KEYWORDS : Bayesian inference, Monte Carlo Markov Chain, Gibbs sampler, Metropolis–Hastingsalgorithm, nonlinear state-space modelsMOTS-CLÉS : Inférence bayésienne, Monte Carlo Chaîne de Markov, échantillonneur de Gibbs, al-gorithme de Metropolis–Hastings, modèles non linéaire à espace d’étatVolume 1 – 2003, pages 1 à 7 – Revue2 Revue – Volume 1 – 2003The purpose of fishery stock assessment is to evaluate the status of a fish stock andto predict how the stock will respond to various exploitation or harvest scenarios. Fishabundance is expressed as the total fish weight (or “biomass”). Increases in the biomassare determined by body growth of individual fish in the population and the addition (or“recruitment”) of the new generations of young fishes (or “recruits”). Those gains mustthen be balanced against the proportion of the population removed by harvesting (calledfishingmortality, F) and other losses due to predation, starvation, or disease (called naturalmortality, M).For each year t = 1 · · ·T of a period of T years, the basic data for stock assessmentsare Ct, the amount of fish caught, and It a measure of the abundance of the stock vulner-able to fishing during this year. The index It is a relative measure of the abundance ofthe stock, here we use the catch per unit of effort (CPUE) (i.e. ratio of fish caught to thetime spent fishing). Using C1:T and I1:T we want to fit a model on these data in order toestimate the biomass Bt for all t, but also to evaluate short–term projections for the stockunder a number of different fishing and management scenarios.Following the work of Meyer and Millar [6], we propose to adopt a state–space rep-resentation of a delay difference model of the biomass dynamic. The delay differencemodel was originally introduced by Deriso [2] and generalized by Schnute [7] (for a gen-eral presentation see Hilborn and Walters [4]). The indexes I1:T will be considered as theobservation process and C1:T as given parameters.One adjusts the parameters of the model to the dataC1:T and I1:T by a MCMC (MonteCarlo–Markov Chain) method, namely a Metropolis–Hastings within Gibbs algorithm.For a general presentation of MCMCmethods one can consult Gilks et al [3], for a surveyof MCMC methods for non linear/non–Gaussian state–space modeling one can consultTanizaki [8].This paper is organized as follows. In the fist section we describe the delay differencestate–space model. The second section introduces the proposed MCMC method : theMetropolis–Hasting method within Gibbs sampling.The results are presented in the lastsection.1. A state–space model for fishery stock assessmentWe are interested in modeling the evolution of the total mass of a fish populationvulnerable to fishery (the biomass) along a given series t = 1 · · ·T of years. At thebeginning of year t, the biomass Bt is :Bt!=!a"kwa Na,t (1)RevueMCMC for state–space models 3where Na,t is the population size at age a in year t, wa is the average weight at age a,and k is the age of recruitment (i.e. the age from which an individual is entered in thebiomass). The first term Rt = wk Nk,t in (1) represents the recruitment that is the part ofthe biomass that has just become susceptible to the fishery.We assume that all recruited fishes are equally vulnerable to natural mortality and tothe fishery. The annual evolution of each cohort is described by :Na,t = st#1 Na#1,t#1 , (2)st#1 denotes the total survival probability in year t ! 1. The evolution of wa is given bythe following weight–at–age growth curve wa = w$ (1 ! e#! (a#a0)) for a " a0, wherea0 is the age for which wa0 = 0, and ! > 0. This gives :wa = (1 + ")wa#1 ! "wa#2 (3)where " = e#! # (0, 1) is the growth rate. One can easily check that (2) and (3) in (1)gives the recurrence :Bt = (1 + ") st#1 Bt#1 ! " st#1 st#2 Bt#2 + Rt ! " st#1 w Rt#1 (4)with w = wk#1/wk # (0, 1). Under the assumption of independence of fishing andnatural mortality, we have st = sMt sFt where sMt is the survival rate to natural mortality,it is supposed to be a constant sM = e#M with M > 0, and sFt is the survival rate tofishery : sFt = Bt#CtBt where Ct is the catch during year t. Then we suppose that therecruitment is constant, i.e. R = Rt. Equation (4) becomes :Bt = (1 + ") e#M Bt!1#Ct!1Bt!1 Bt#1 ! " e#2M Bt!1#Ct!1Bt!1Bt!2#Ct!2Bt!2Bt#2+ R"1 ! " e#M w Bt!1#Ct!1Bt!1#.We suppose that the catch is null for t < 0, and that Bt is at equilibrium before timet = 1, i.e. B1 = K (the virgin biomass), and :B2 = (1 + "! " e#M ) e#M (B1 ! C1) + R"1 ! " e#M w B1#C1B1#.The (deterministic) observation equation is It = q Bt for t = 1 · · ·T , where It is arelative biomass index and q is the “catchability” coefficient.Revue4 Revue – Volume 1 – 2003State and observation error processes are modelized as a multiplicative lognormal i.i.d.processes (Meyer and Millar [6] propose a conditional additive Gaussian noises). We endwith the following model :B1 = K $ e%"2W W1 , (5)B2 = F2(R, B1) $ e%"2W W2 , (6)Bt = Ft(R, Bt#1, Bt#2) $ e%"2W Wt , t = 3 · · ·T , (7)It =$q Bt%$ e%"2V Vt , t = 1 · · ·T , (8)W1:T and V1:T are N(0, 1) Gaussian noises, andF2(R,B1) = (1 + ! ! ! e!M ) e!M (B1 ! C1) + R 1 ! ! e!M w B1!C1B1 ,Ft(R, Bt!1, Bt!2) = (1 + !) e!M Bt!1!Ct!1Bt!1Bt!1 ! ! e!2M Bt!1!Ct!1Bt!1Bt!2!Ct!2Bt!2Bt!2+ R 1 ! ! e!M w Bt!1!Ct!1Bt!1 .In this application the parameters ", M , w are supposed to be given. The parametersK , R, q, #2W , #2V are supposed to be random variables. These parameters, and the noiseprocessesW1:T and V1:T are supposed to be mutually independent. The a priori laws ofK , R, q, #2W , #2V are specified later.This model is a good tradeoff between surplus production models and age–structuredmodels : from one hand surplus production models are too simplistic (in these modelsall the components of the evolution of the biomass are aggregated into a unique one),and from the over hand age–structured models are more sophisticated but rely on age–structured catch data which are not available is most applications.2. MCMC methodThe catches Ct are given. The parameters ", $, and M are supposed to be given. Inthis first work, we also suppose that the variance parameters #2W , #2V are known, even ifit is a restrictive hypothesis. Hence we want to inferK , R, q, B1:T from the observationsI1:T . For the sake of notational simplicity let :X1:P = [K, R, q, B1:T ] with P = T + 3 .Our goal is to sample from the a posteriori target density1 :p(X1:P |I1:T ) . (9)1. For convenience, X will denote a random variable and its realizations, p(X) will denote the densityof this rv (it will also correspond to either a density or a probability measure), p(X|Y ) will denote theRevueMCMC for state–space models 5Because sampling from this joint conditional density is not straightforward we propose touse a MCMC algorithm, i.e. to build a Markov chain {X (k)1:P}k"0 whose ergodic invariantmeasure admits the target density (9) so that simulating this chain will give samples fromthe target density.Require: for each component p : the proposal sampler "p( · |X\p, I1:T ) and the likelihood#p(Xp, X\p, I1:T ) (cf. (11)).Ensure: sampling from the target density p(X1:P |I1:T ) (cf. (9))k " 0choose initial guess X1:Ploopfor p = 1, . . . , P doX "p " "p( · |X\p) {sample the new candidate}$ = #p(X"p|X\p, I1:T )/#p(Xp|X\p, I1:T ) {cf. (12)}U # U [0, 1]if U $ $ thenXp " X "pend ifend fork " k + 1if k % kburn–in thenreturnX1:P {after the burn–in period}end ifend loopAlgorithm 1: Metropolis–Hastings within Gibbs samplerWe apply the Gibbs sampler : we iteratively sample from each of the marginal condi-tional densities ofXp given other componentsX{1,...,P}\{p} (and I1:T ) i.e. from :p(Xp|X\p, I1:T ) with \ p!= {1, . . . , P} \ {p} (10)for p = 1, . . . , P . Even if these densities admit explicit analytical representations, wecannot sample directly from them. For this reason, we use a Metropolis–Hasting algo-rithm within the Gibbs sampling loop. Each of the previous densities will be decomposedas the product of a density we know how to sample from, and of a likelihood that we cancompute explicitly (up to a multiplicative constant) :p(Xp|X\p, I1:T ) & %p(Xp|X\p, I1:T )& '( )proposal density$&p(Xp, X\p, I1:T )& '( )likelihood. (11)conditional density law of X given Y (we should write x !" pX|Y (x|y) for the conditional density of Xgiven Y = y). The meaning will be clear from the context.Revue6 Revue – Volume 1 – 2003Starting from a given configurationX (k)1:P we sample the next configurationX(k+1)1:P asfollow : let X(k+1)1:P ' X(k)1:P then for each component p # {1 : P} we sample a new pthconfigurationX %p ( %p( · |X(k+1)\p , I1:T ) ,then with probability 1 ) ' where :'!=#p(X"p,X(k+1)\p ,I1:T )#p(X(k)p ,X(k+1)\p ,I1:T ), (12)we accept X %p, i.e. X(k+1)p = X %p (with probability 1 ! 1 ) ' the pth component isnot changed). All the components are updated in the same way. After a certain period(burn-in) the Gibbs sampler reaches its stationary behavior and produce samples from theposterior density. The detailed procedure is summarized in Algorithm 1. The marginalconditional densities and their decomposition between proposal densities and likelihoodfunctions are partly detailed in the appendix.3. Results and commentsWe performed this algorithm on simulated data which match the data used in Meyerand Millar [6] (themselves taken from Kimura et al [5]). The stationarity of the sampleris not reached on the component q in less than millions of iterations, so we chose to frozethis parameter to its true value (i.e. supposed that it is known). Then 250000 iterations ofthe Gibbs sampler were performed (with a burn-in period of 1000 iterations).For the simulation and the MCMC algorithm we used the following values : " =1, $ = 0, M = 0.6, #2W = #2V = 0.01 (the parameter q is also supposed to beknown q = 9.39). Moreover for the simulation we use the values R = 218.5, andK = R 1#$ e!M %1#(1+$) e!M +$ e!2M * 1073.3 (which corresponds to the virgin biomass, i.e. theasymptotic solution of (4) without catch), the biomass time series Bt=1934···1967 and theCPUE one It=1934···1967 are simulated through the system (5–8).As illustrated by Fig. 1, the match between the parameters used in the simulation andthe a posteriori law inferred by the Gibbs sampler (cf. Figure 1) is fairly good. Thecorrespondence between the “true” (simulated) biomass time series and the correspondinga posteriori law derived from the Gibbs sampler is also good.The main drawback of this approach is that it is computationally too slow and cum-bersome, so it seams necessary to use more elaboratedMCMC techniques like interactingMetropolis procedures (Del Moral and Doucet [1]). Moreover, it is necessary to com-pare this approach to a non Bayesian one (e.g. to approximate the maximum likelihoodestimator with an EM algorithm).RevueMCMC for state–space models 70.5 1 1.5 2x 10560080010001200140016001800MCMC iterations K600 800 1000 1200 1400 1600 1800empirical a posteriori law true value0.5 1 1.5 2x 105200220240260280MCMC iterations R200 220 240 260 280empirical a posteriori law true value0.5 1 1.5 2x 105500600700800MCMC iterations B1958500 600 700 800empirical a posteriori law true valueyear  tmillions of pounds1935 1940 1945 1950 1955 1960 1965500100015002000 BtFigure 1. For the components K, R, B1958, the iterations of the MCMC procedure aredisplayed on the left and the resulting empirical a posteriori laws (together with the truevalues of the parameters) on the right. Below is displayed the biomass time seriesBt with,for each year t, the corresponding empirical a posteriori law represented in grey levels.AcknowledgementsM. Rakotozafy’s work has been partially supported by the French Ministry of ForeignAffairs (French embassy in Madagascar), part of this work was carried out during a stayat IRISA.4. ReferencesP. Del Moral and A. Doucet. On a class of genealogical and interacting metropolis models. In J.-M. Morel, F. Takens, and B. Teissier, editors, Séminaire de Probabilités XXXVII, volume 1832of Lecture Notes in Mathematics, pages 415 – 446. Springer–Verlag, 2003.Revue8 Revue – Volume 1 – 2003R. B. Deriso. Harvesting strategies and parameter estimation for an age-structured model. Cana-dian Journal of Fisheries and Aquatic Sciences, 37:268–282, 1980.W. Gilks, S. Richardson, and D. Spiegelhalter, editors. Markov Chain Monte Carlo in Practice.Interdisciplinary Statistics. Chapman and Hall, 1996.R. Hilborn and C. Walters. Quantitative Fisheries Stock Assessment : Choice, Dynamics andUncertainty. Chapman and Hall, 1992.D.K. Kimura, J.W. Balsiger, and D.H. Ito. Kalman filtering the delay difference equation: practicalapproaches and simulations. Fishery Bulletin, 94(4):678–691, 1996.R. Meyer and R.B. Millar. Bayesian stock assessment using a state-space implementation of thedelay difference model. Canadian Journal of Fisheries and Aquatic Sciences, 56:37–52, 1999.J. Schnute. A general theory for the analysis of catch and effort data. Canadian Journal ofFisheries and Aquatic Sciences, 42:414–429, 1985.H. Tanizaki. Non linear and non-Gaussian state–space modeling with Monte Carlo techniques : Asurvey and comparative study. In C.R. Rao and D.N. Shanbhag, editors, Handbook of Statistics21: Stochastic Processes: Modeling and Simulation, pages 871–929. North–Holland, 2003.Appendix : the marginal conditional densitiesThe marginal conditional densities (10) are :➀ p(K|R, q, B1:T , I1:T )➁ p(R|K, q, B1:T , I1:T )➂ p(q|K, R, B1:T , I1:T )➃t p(Bt|K, R, q, B\t, I1:T ) , t = 1 · · ·TFor each of these marginal densities we specify the decomposition (11) between the pro-posal density and the likelihood function :➀ & %(K) ' p(B1|K)➁ & %(R) ' p(B2|R, B1) Tt=3 p(Bt|R, Bt!1, Bt!2)➂ & %(q) ' Tt=1 p(It|q, Bt)➃1 & p(B1|K) ' p(B2|R, B1) p(B3|R, B2, B1) p(I1|q, B1)➃2 & p(B2|R, B1) ' p(B3|R, B2, B1) p(B4|R, B3, B2) p(I2|q, B2)➃t & p(Bt|R, Bt!1, Bt!2) ' p(Bt+1|R, Bt, Bt!1) p(Bt+2|R, Bt+1, Bt)p(It|q, Bt)➃T!1 & p(BT!1|R, BT!2, BT!3) ' p(BT |R, BT!1, BT!2) p(IT!1|q, BT!1)➃T & p(BT |R, BT!1, BT!2) ' p(IT |q, BT )proposal density likelihoodand the prior densities on the parameters (whose choice will be detailed elsewhere) are :((K) ( U [400, 4000], ((R) ( U [0, 2000], ((q) ( U [0, 100] (uniform laws).Revue

{MCMC} for non linear/non {Gaussian} state-space models: Application to fishery stock assessment

https://hal.inria.fr/hal-00652092/document

{MCMC} for non linear/non {Gaussian} state-space models: Application to fishery stock assessment

Abstract

Similar works

Full text

Available Versions

INRIA a CCSD electronic archive server