On the stochastic approach to marine population dynamics

The purpose of this article is to deepen and structure the statistical basis of marine population dynamics. The starting point is the correspondence between the concepts of mortality, survival and lifetime distribution. This is the kernel of the possibilities that survival analysis techniques offer to marine population dynamics. A rigorous definition of survival and mortality based on their properties and their probabilistic versions is briefly presented. Some well established models for lifetime distribution, which generalise the usual simple exponential distribution, might be used with their corresponding survivals and mortalities. A critical review of some published models is also made, including original models proposed in the way opened by Caddy (1991) and Sparholt (1990), which allow for a continuously decreasing natural mortality. Considering these elements, the pure death process dealt with in the literature is used as a theoretical basis for the evolution of a marine cohort. The elaboration of this process is based on Chiang ́s study of the probability distribution of the life table (Chiang, 1960) and provides specific structured models for stock evolution as a Markovian process. These models may introduce new ideas in the line of thinking developed by Gudmundsson (1987) and Sampson (1990) in order to model the evolution of a marine cohort by stochastic processes. The suitable approximation of these processes by means of Gaussian processes may allow theoretical and computational multivariate Gaussian analysis to be applied to the probabilistic treatment of fisheries issues. As a consequence, the necessary catch equation appears as a stochastic integral with respect to the mentioned Markovian process of the stock. The solution of this equation is available when the mortalities are proportional, hence the use of the proportional hazards model (Cox, 1959). The assumption of these proportional mortalities leads naturally to the construction of a survival model based on the Weibull distribution for the population lifetime. Finally, the Weibull survival model is elaborated in order to obtain some reference parameters that are useful for management purposes. This section does not deal exhaustively with the biological and fishery reference parameters covered in the specialised monographs (Caddy and Mahon, 1996; Cadima, 2000). We focused our work in two directions. Firstly, the principal tools generating the usual reference parameters were adapted to the proposed Weibull model. This is the case of biomass per recruit and yield per recruit, which generate some of the important reference points used for management purposes, such as the FMSY, F0.1, Fmed. They also provide important and useful concepts such as virgin biomass and overexploitation growth. For this adaptation, it was necessary to previously adapt the critical age as well as the overall natural, fishing and total mortality rates. Secondly, we analysed some indices broadly used in all population dynamics (including human populations) but only marginally dealt with in fishery science, such as life expectancy, mean residual lifetime and median survival time. These parameters are redundant with mortality rates in the classical exponential model, but are not so trivial in a more general framework.


INTRODUCTION
Fisheries science has been developed through deterministic models with a simple mathematical background that can be understood by fishery biologists.There exist methods that provide appropriate computer tools for diagnosing the state of exploitation of resources, mainly by comparing the current situation with the maximum sustainable yield.The works by Gulland (1974Gulland ( , 1977Gulland ( , 1983)), Ricker (1975) and Beverton and Holt (1957) are crucial steps in this initial deterministic approach.In these basic studies the authors propose constant mortality rates that are independent of age, or step functions over unit intervals of age that are coherent with the exponential evolution of the size of a cohort.One of the most popular evaluation procedures, virtual population analysis (VPA) and all its derivations, use this axiom in order to establish the assessment through successive retrospective evaluations of the stock.The assessment is a process of inference based on the mortality models and relating two demographic structures: that of the observed catches of the cohort (or pseudo-cohort) and that of the unknown population.
The analysis of the two components of the total mortality rate, natural and fishing mortality, plays an essential role in the elaboration of this idea.The difficulties of exploring the sea in a virgin state imply that, in most cases, natural mortality is not experimentally assessed but indirectly evaluated and assumed to be constant in the classical deterministic approach.Caddy (1991), Sparholt (1990) and Chen and Watanabe (1989) have proposed alternatives to this assumption of constant natural mortality and, in general, as in Abella et al. (1997), these models have been used to build a vector of natural mortalities (or natural mortality as a step function) in order to improve the assessment based on deterministic methods.In consequence, very little work has been done relating the fundamental models of marine population dynamics to the underlying lifetime distribution because, according to the previous description, it has been systematically assumed to be piecewise exponential corresponding to a mortality rate formalised as a step function.Although the essential random character of the phenomena and of the available data has not been considered in the classical deterministic approach, this scenario is recently changing, as Buckland et al. (2000) pointed out.Already in the works by Gudmundsson (1987) and Sampson (1990), the evolution of the stock of a cohort is formalised as a stochastic process.
In the present work, the stochastic nature of marine population dynamics is structured on the basis of the possible lifetime distributions and their corresponding mortality and survival models.The survival models form the basis for considering the evolution of the stock of the cohort and its cumulative catches as stochastic Markovian processes in the way established by Chiang (1960) for the general life table of any population.On one hand, this structure provides new stock assessment methods that improve the classical ones.On the other hand, it establishes a promising bridge between marine population dynamics and three important and powerful stochastic fields and associated tools: stochastic processes, survival analysis and multivariate analysis.SCI. MAR.,71(1), March 2007, 153-174. ISSN: 0214-8358 154 • E. FERRANDIS de estos procesos a través de procesos gaussianos, permite la aplicación del importante cuerpo teórico y computacional del análisis gaussiano multivariante para el tratamiento probabilístico de problemas relevantes relacionados con la actividad pesquera.Como una consecuencia, la necesaria "ecuación de captura" aparece como una integral estocástica con respecto al mencionado proceso markoviano del stock.Y la solución de dicha ecuación es inmediata cuando las mortalidades asociadas a las distintas causas de muerte son proporcionales.De aquí se deduce la importancia y posibilidades del conjunto de técnicas conocido como "riesgos proporcionales" y establecido fundamentalmente por Cox (1959).Asumir estas mortalidades proporcionales, conduce de forma natural a considerar un modelo de supervivencia basado en la distribución de Weibull para el tiempo de vida de la población.Finalmente, se elabora dicho modelo de supervivencia de Weibull para obtener algunos parámetros de referencia utilizados en la gestión de los recursos marinos explotados.Esta sección no trata de abarcar exhaustivamente los parámetros de referencia biológico-pesqueros que se tratan en los específicos textos de referencia (Caddy y Mahon, 1996;Cadima, 2000) y se focaliza en dos direcciones.En primer lugar los principales parámetros usualmente utilizados en la gestión de los recursos: es el caso de la biomasa por recluta y el rendimiento por recluta que generan algunos de los importantes puntos de referencia utilizados en la gestión, tales como los F MSY , F 0.1 , F med , así como importantes conceptos como los de biomasa virgen o sobreexplotación de crecimiento.Para esta generalización, ha sido necesaria la adaptación previa de la edad crítica, así como las de tasas globales de mortalidad.En segundo lugar, se analizan algunos índices de uso extensivo en dinámica de poblaciones (incluida la humana) pero tratadas de forma marginal por la ciencia pesquera: es el caso de la esperanza de vida, la vida mediana y el tiempo de vida residual.Estos parámetros son redundantes con las tasas de mortalidad en el modelo exponencial clásico pero no resultan inmediatos en un modelo más general, como el propuesto.

THE FUNDAMENTAL TRIANGLE
Let T be the lifetime (duration of life) of the population, considered as a random, real valued, nonnegative and absolutely continuous variable.Let f(t) and F(t) be its density and distribution functions, respectively.Then, the survival function is defined by (1.1) or (1.2) and the total mortality rate or "force of mortality", is expressed as (1.3) which is the "hazard function" or "hazard rate" in classical survival analysis (Smith, 2002).Therefore, Z(t) is the instantaneous rate of death at age t, given that the individual survives up to age t.Or, in the conceptual framework systematised by Cadima (2000), the total mortality rate should be the instantaneous relative rate of the survival (with a positive sign).
The integral of the mortality (1.4) is the "cumulative hazard" (Altman, 1999) or the "integrated hazard" (Dobson, 2002).Its relationship with the survival function, given by the integration of Equation (1.3), becomes (1.5) By the well known properties of distribution and survival functions: Then we have (1.6)Let t Max denote the limit of the attainable ages of the population (it can be finite or infinite, known or unknown).Then (1.7) (1.8) Therefore, the effective range of the ages will be [0, t Max ).However, this is equivalent to considering the infinite interval [0, ∞) as the range of the ages and a null mortality rate and survival, S(t) = Z(t) = 0, for every age greater than t Max .
The integrability conditions of the mortality rate imply (Munroe, 1953: 191) that (1.9) which, in turn, implies the continuity of the survival function S(t).
Also, the expression (1.5) of the survival as the exponent of an integral implies (Munroe, 1953: 268;Dieudonné, 1963: 159) that the survival must have a derivative almost everywhere (a.e.).This means that the derivative may not exist at a set of measure zero, like a set of isolated points.In fact, this is the usual case considered in classical marine population dynamics, when the mortality rate is a step function and the survival is not differentiable at the discontinuities of the mortality rate.
We can now summarise the properties of mortality and survival.
The mortality rate must be non-negative, real valued, integrable (a member of the space L 1 of Lebesgue-integrable functions) in every finite interval [0,t] bounded by an age t lower than the maximum attainable age t Max , but non-integrable in any interval containing the [0, t Max ) interval, in particular the interval [0,∞): The first condition is related to the decreasing (at least non-increasing) behaviour of the survival.The second condition guarantees the existence of the survival and its initial unit value.The third condition implies that, as the age increases, the limit of the survival is null.Any function verifying properties (1.10) is an admissible mortality rate whose survival function is given by Equation (1.5).
The survival function is real valued, bounded by the [0,1] interval, non-negative and continuous, defined on all the positive ages, the [0,∞) interval.The space of functions with these properties is usually designated by C [0,1] [0,∞) (Dieudonné, 1963: Chapter VII).In addition it is non-increasing, with an initial value of 1 and tending to 0 as the age increases.Finally, it is differentiable almost everywhere.That is: (1.11) Any function verifying these properties is an admissible survival function.

THE ANALYSIS OF THE MORTALITY RATE: INDEPENDENT AND COMPETITIVE RISKS
The mortality rate is generally expressed as a sum of natural (M) and fishing (F) components: (2.1) To say that the natural mortality rate, M(t), should be the total mortality rate in absence of exploitation, is equivalent to assuming that the two causes of death act independently.
In general, if there are k risks of death R i (i = 1,2,…,k) which act simultaneously on each individ-ual in a population, Z i (t), S i (t) are the corresponding mortality rates (or risk functions) and survival of the lifetime which would be applied if R i were the only risk present, and Z(t), S(t) are the total mortality rate and survival, then the three following conditions are equivalent (Cox, 1959, andDavid, 1970): This hypothesis of the independence between natural and fishing mortality and the corresponding analysis of the mortality rate (2.1) is usually assumed in fisheries science, although it may be far from reality, particularly for a multispecies fishery, and even more particularly for demersal resources.The fishery gears could represent a certain "competition" for the natural predators and a deflation factor for their populations, because they will often be captured too.Consequently, the dependence between fishing and natural mortality should be quite realistic.In this case, these two causes of death should be considered as dependent or "competitive" risks (Chiang, 1970(Chiang, , 1991;;David, 1970;Gail, 1975).

EXAMPLE MODELS
In the usual applications of virtual population analysis, the mortality rates are considered as step functions.The age intervals defining them are in most cases of one year, with the exception of the last.In these cases, the non-integrability of the mortality rate in the whole interval [0,∞) implies that the last age interval must be a plus-class and in this last interval the total mortality rate must be strictly positive.Otherwise, the proposed step function must be integrable in [0,∞) and could not be considered as an admissible total mortality rate.As a particular case, a constant positive function on the whole interval [0,∞), as it is frequently considered the natural mortality rate, is an admissible mortality model that verifies conditions (1.10).
As an example we shall consider the "reciprocal function" proposed by Caddy (1991) ´( ) . .[ , ) The integral of this function is The function M(t) is not integrable in any finite or infinite interval including the initial age of the cohort (the integral has an infinite limit when the age approaches zero).Conditions (1.10) are not satisfied, so it cannot be an admissible model for a natural mortality rate.
An adequate and admissible alternative to Equation (3.1) would be: Now the integral is Hence, all conditions (1.10) are verified and Equation (3.2) defines an admissible mortality rate.
Another admissible mortality decreasing with age, and therefore a possible model for the natural mortality rate, could be: with a and b > 0 whose integral is which also verifies all conditions (1.10).Therefore, Equation (3.3) defines an admissible mortality rate whose corresponding survival function is: Caddy (1991) refers to an unpublished suggestion made by R.J.H. Beverton in 1991 for a derivation of the reciprocal model leading to the following expression: (3.4) with a, b and r > 0, interpreting r as an initial considered age.This suggestion by Beverton is indeed an extension of the reciprocal model and is an admissible mortality model whose integral is and which verifies conditions (1.10).
Another natural mortality model cited by Caddy (1991) and due to Sparholt (1990) is (3.5) which is also a decreasing function of age, converging to the asymptotic null value as the age increases.
The integral of this function is which is integrable on the whole range of ages, i.e. the [0,∞) interval.Therefore, the function (3.5) cannot be by itself an admissible mortality model but only a component of a more complex model.

THE ELABORATION OF THE BEVERTON AND HOLT AND SPARHOLT MODELS:
PROPORTIONAL HAZARDS Beverton and Holt (1957) proposed a very simple and yet powerful model for the natural mortality rate, as the basic hypothesis for the stock-recruitment relationship.The mortality rate, Z(t), should have two components, one constant and the other depending on the number of surviving individuals.We adapt this idea considering the proportion of survivals, i.e. the survival function S(t).The model has an initial value Z 0 = Z(0) and converges asymptotically to the adult mortality Z a = lim t→∞ Z(t).The acceleration of this convergence is dominated by the shape parameter α.For a positive shape, α>0 (4.1) defines a mortality decreasing with age as does the survival S(t).
When the shape parameter is zero, we have the constant mortality rate and the simple exponential distribution for the lifetime of the population.
Combining this expression with Expression (1.3) of total mortality, we have the following differential equation which may be expressed as

S t bt
If the adult mortality, Z a , is null, the integration of (4.2) gives (4.3) which is a particular case of the admissible model (3.3).
Otherwise, if Z a >0, we can rescale the shape parameter by expressing the power of the survival as α/Z a .Then, the corresponding relation between mortality and survival becomes (4.4)Now, the differential equation ( 4  In general, the mortality model (4.5) is compatible with proportional mortality rates corresponding to different causes of death acting independently.Indeed, if a partial component of the mortality (i.e. the natural mortality, M(t)) presents a mortality rate of the same type (with a shape parameter β) and proportional to the total mortality, Hence and finally Thus, the proportional mortalities imply a common shape parameter and the total mortality may split as the sum of proportional fishing and natural mortalities with the same shape parameter and acting independently.This is formalised in terms of mortality (4.11) or in terms of survival (4.12) In the same way, an admissible extension of the Sparholt model for natural mortality is where, as in the previous model, M 0 = M(0), is the initial mortality, M a = lim t→0 M(t) is the asymptotic or adult mortality and α is a shape parameter which controls the more or less accelerated transition from the initial to the adult mortality or a measure of the premature character of the mortality.This model is also coherent with the possible proportionality of the mortalities due to the different causes of death, if we consider the fishing and total mortalities following the same mortality pattern, with the same shape parameter and express the proportional factor in terms of an assumed constant exploitation rate, E.  (4.16)This proportionality between fishing, natural and total mortalities agrees with the idea of considering the fishery as another element of the ecosystem, a powerful predator whose effects should be coherent and related (proportional) to the other causes of mortality.This hypothesis, as will be seen, is essential in order to obtain an available catch equation.It may be assumed in the present generalised models (4.4) and (4.13), as well as in the Weibull model proposed below.In all these models, the natural mortality rate is assumed to be the same in the absence or in the presence of the fishery and both causes of mortality, natural and fishing, are considered to act independently-an independence hypothesis that is a common assumption in fishery science.
If an age of first capture, t c , is considered, a sim-ple and canonical model leads to the following battery of lifetime distributions corresponding to the unexploited phase, subject only to natural mortality, and the exploited phase, subject to total mortality.Hence (4.17) where This decomposition of the range of ages (into two or more subintervals) is particularly useful when the conditional survival (4.18) is available and the mortality rates are proportional in each interval, as is the case with the elaborated Beverton and Holt and Sparholt models.

THE LIFETIME, SURVIVAL AND EVOLUTION OF THE STOCK
We assume a marine cohort with an initial population size N 0 in its initial age t = 0, a function Z(t) for the mortality rate and a survival function S(t).Let T be the lifetime of the population and T i ,i = 1,2,...,N 0 the lifetimes of the different individuals of the cohort.
The size of the stock at age t, N(t), is the number of survivors at that age and may be expressed as (5.1) where I [t,∞] is the indicator function of the interval [t,∞], defined, for any set A, by If we assume the independence of the lifetimes of the individuals of the cohort, N(t) has the binomial distribution with parameters N 0 for the sample size and Pr (T ≥t) = S(t) for the probability of success.Its expectation, m(t), and variance, v(t), are (5.2) As a consequence of (1.3), we have For any two ages, t 1 ≤t 2 , the covariance between the cohort sizes may be obtained as follows and hence (5.4) The mean and variance of the number of deaths in any interval of ages [t 1 ,t 2 ], N(t 1 ) -N(t 2 ) are: (5.5)An alternative definition of the process may be based on the fact that, for a cohort of initial population size N 0 , the expected population decrease at age t, i.e. the expected number of deaths is N 0 (1-S(t)).If we consider the number of deaths, N 0 -N(t), as a random variable with a Poisson distribution, its expectation, variance and parameter, v*(t) = N 0 (1-S(t)), coincide and reproduce the variance of the size of the cohort N(t): var(N(t)) = var(N 0 -N(t)).
If the numbers of deaths in disjoint intervals is considered independent, the number of deaths of the cohort at the age is a Poisson process (non-stationary).
With this definition, the mean, variance and covari-ance functions of the size of the cohort at age t are: (5.6) In this Poisson process, the mean and variance of the number of deaths, N(t 1 ) -N(t 2 ), in any interval of ages [t 1 ,t 2 ] are: (5.7)
In the case of the binomial process, due to its construction the conditional distribution of N(t), given {N(t 1 ),N(t 2 ),...,N(t k )} has the binomial distribution or which verifies the Markovian property: "the future depends on the past only through the present".
And this result is coherent since Using the binomial probability function it is easily verified that The mean and variance of the conditional distributions are: ´( ) which gives (6.1) and similarly (6.2) Notice that Equation (6.1) is the so-called "stock equation" or "fundamental equation of population dynamics".In the case of a constant mortality rate Z in the interval considered, we have (6.3)A necessary and sufficient condition for the conditional distribution of the process reproduces the fundamental equation of population dynamics as a conditional expectation, is that for every t i < t j we have covariance = variance (6.4) Indeed: c being the regression coefficient, Taking into account Equation (5.2), we have Hence, the condition that becomes necessary and sufficient.
The evolution of the size of the cohort, N(t), with its mean, variance and covariance functions (5.2) and (5.4) is a binomial and Markovian process whose conditional distributions reproduce the fundamental equation of population dynamics as a conditional expectation.
If the number of deaths is considered as a nonstationary Poisson process, condition (6.4) is not satisfied and the stock equation only holds for the expected stock size of the cohort: (6.5) In this case, the conditional variance is equivalent to the variance of the increments, as a consequence of the independence of the increments: The independence of the increments of a Poisson process implies that here again the Markovian property is verified.

THE LOGNORMAL VERSION OF THE STOCK PROCESS
In a deterministic approach, the stock of a cohort of marine fishes under exploitation, N(t), submitted to a total mortality rate Z(t) and starting from an initial population size N 0, evolves over time (or age) t according to the differential equation A way to model this evolution through a stochastic process is presented by Oksendal (1985).The mortality rate is often randomised by introducing a "noise" where the "noise" is formally proportional to a white noise, X(t) (normal and independently distributed): As an example, the simulation of stocks made by the ICES Methods Group (ICES, 1988) takes a constant natural mortality and obtains the randomisation of the fishing mortality rates, F, by a lognormal var( ( ) where λ is a constant, X the standard normal distribution, F the "expected" fishing mortality rate and F* the achieved or "realised" one.The error, F*-F, is called the "process error" and it is attributed to climate or other external factors that may alter a previously decided fishing effort.This randomisation is equivalent to and taking the first approximation of the Taylor series development of the exponential function, we get This randomisation of the mortality rate originates the stochastic differential equation W(t) being the standard Brownian motion and b being considered as a constant in the simplest modelling.The solution of this stochastic equation is (Oksendal, 1985) (7.3) In this case, it is the natural logarithm of the stock that has a normal distribution and is expressed as a generalised Brownian motion.The independence of the increments of the Brownian motion guarantees the Markovian character of the lognormal stock process.
This lognormal version is a more complex process than the previous binomial and Poisson one, as it incorporates the noise parameter(s) "b", and it has the following expectations and covariances: (7.4) The covariance function verifies condition (6.4).Hence this model provides the stock equation (6.1) as a conditional expectation.

THE CATCH PROCESS
Let us assume the usual decomposition (2.1) of the total mortality as a sum of two independent components: natural (M) and fishing (F) mortalities.
The simplest expression of the catch process in a cohort, C(t) or C[0,t], as the random number of individuals of up to age t caught by the fishery, is through the stochastic integral (8.1) or, for any interval of ages [t 1 ,t 2 ] (8.2) Considering the binomial process defined before, the conditional expectations of the increments in an age interval, [h,t], are bounded by the following monotonous function and The function 1 -S(t)/ S(t) is continuous and monotone increasing.Hence, the binomial process verifies the conditions established by McShane (1969McShane ( , 1974)), and are nigh-martingales in the sense defined by Young (1970Young ( , 1974)).These are sufficient conditions for the existence of the stochastic integral (8.1) or (8.2).In fact, the stochastic integral exists in much more general conditions than those considered "noise" " in (8.1) or (8.2), with an integrand which is a real valued and bounded function, as it is the exploitation rate.The existence includes all the stochastic processes whose trajectories are integrable.Hence, it allows the mortality rates to be randomised, as is shown when the lognormal version of the stock process is dealt with.The stochastic integral (8.1) always exists and its solution is the catch equation.
Moreover, by the properties of the integral, for any real valued function, f(t), the expectation of the integral is This means that any concept (for instance, the yield per recruit) derived in the deterministic case from the catch in number becomes the corresponding expected value (the same concept in mean) when the catch in number is considered as a stochastic process.
The catch equation has an immediate solution if the exploitation rate F(t)/Z(t) is constant in the considered interval of ages.It is the same to say that the mortality rates are proportional within this interval.In this case, the solution of the integral is (8.3)Then (8.4)If the stock process verifies the covariance condition (6.4), as the binomial and the lognormal processes do, the catch equation becomes (8.5) In the simplest case, when the mortality rates are constant in the interval of ages considered, the catch equation leads to the well-known expression (8.6) which is again a conditional expectation.This expression is usually simplified in its deterministic version as: (8.7)An analogous development corresponds to the considered Poisson process.

and
The function 1-S(t) is continuous and monotone increasing.Hence, this non-stationary Poisson process verifies the above-mentioned existence conditions for the catch process.In this case, the catch equation appears as an unconditional expectation. (8.7)

GAUSSIAN APPROXIMATIONS
A normal approximation to the binomial distribution is justified by the great expectations present when marine populations are dealt with, leading to a Gaussian process with mean and covariance functions given by (5.2) and (5.4) respectively.
The normal approximation to the lognormal process has the mean and covariance functions given by (7.4).In both cases the mean and covariance structures may be integrated as follows (9.1)where (9.2) 1 for the lognormal process s for the binomial process The function B(t), in both cases, is monotone increasing with B(0) = 0.
These processes verify the Kolmogorov conditions for stochastic processes.These conditions refer to permutations and subsets of the finite dimensional distributions that in the case of multivariate normal, remain multivariate normal.It is therefore only necessary to show the coherence of the covariance structure (9.1), or in mathematical terms, that the covariance matrix corresponding to any finite dimensional distribution is semi-positive definite.For any set of ordered ages t 1 <t 2 <...<t k , the corresponding normal multivariate distribution is with and , for i≤j.
The determinant of the covariance matrix is: The proof of this statement is deduced by induction on k.It is obvious for k = 1.Now, if we subtract the first row of ∑ multiplied by S(t i )/S(t 1 ) from the successive i th row, we reduce the determinant to a k-1 order.
where and the statement is proved by applying the proposed expression to the determinant of ∑ * , which is analogue to ∑.As B(0)<B(t 1 )<...<B(t k ), the determinant remains positive for any diagonal box of the proposed covariance matrix ∑ and is therefore a positive definite matrix.
In practical applications it is decisively important to manipulate the inverse, ∑ -1 , of the covariance matrix.This inverse may be expressed by the tridiagonal matrix: (9.4) where as it is verified that .These expressions of the determinant and the inverse of the covariance matrix enable and simplify the possible applications such as likelihood estimation.
Defined by these mean and covariance structures, for any set of ordered ages t 1 <t 2 <...<t k <t , it is easy to get the expectations and variances of the conditional distributions, N(t) |{N(t 1 ),N(t 2 ),...,N(t k )}: (9.5) (9.6)Which implies that the well-defined models are Markovian processes.
The Gaussian approximation to the alternative nonstationary Poisson process leads to a Gaussian process of independent ("orthogonal") increments ("additive processes"), whose distributions are normal.
var ( ) | ( ), ( ),..., ( ) var ( These processes are known as the generalised Brownian motion, that is, a Brownian motion with non-linear mean and variance functions, a particular case of Gaussian martingales (Doob, 1953;Yeh, 1973).They are a natural extension of the work by Einstein (1905) and Wiener (1923) and constitute the most general family of continuous additive processes.
The additive character of the process, i.e. the above-mentioned independence of the increments, constitutes the weak point of the generalised Brownian motion but is a practical advantage when the parameters of the model are estimated by maximum likelihood.
In the present case, the increments of the process are decrements in the population of the considered cohort, i.e. the number of dead individuals.
Note that from (4.7) Thus, this process may be improved with a greater versatility in allowing a more flexible relation between means and variances by introducing the concept of an "aggregation function", in the sense given to this term by Pielou (1969).
This aggregation function, a(t), may be defined as the ratio between variances and means of the increments of the process (9.8) Starting from the "mean generating function", m(t) = N 0 S(t), a model for the aggregation function produces the corresponding "variance generating function" and the corresponding version of the generalised Brownian motion defined by the finite dimensional distributions where the coordinates of the mean vector and the elements of the covariance matrix are given by the expressions and 10.GENERAL PROPORTIONAL HAZARD MODELS David (1970), elaborating sufficient conditions for the proportionality between mortality rates, proposes three types of two-parameter lifetime distributions that should be coherent with this property: Among these distributions, only the Weibull model generalises the usual exponential model broadly applied in marine population dynamics.Furthermore, the Weibull model can be used for a mortality rate with high values for initial ages and continuously decreasing values for increasing age.The Weibull distribution for the lifetime is thus a "natural" theoretical scenario.In fact, Smith (2002) shows that it represents a very general model coherent with the proportional hazards (between total and natural mortality) implied by the constant exploitation rate.

THE WEIBULL DISTRIBUTION AND THE EXTENSION OF CLASSICAL HYPOTHESES
The Weibull survival (with its corresponding lifetime distribution) is given by (11.1)where Z and α are positive real values.
The usual exponential distribution corresponds to the particular case in which the exponent α is one.
The exponent α determines the "shape" of the survival (mainly in the early ages), the density functions (with a finite maximum for α >1) and the mortality rate (decreasing with age if α < 1, constant for

S t T t e Zt
the exponential case and increasing with age if α >1).The coefficient Z expresses the relative intensity of the mortality, and consequently the relative diminution of the survival.
It is appropriate to identify Z as a mortality coefficient, which is not to be confused with a constant mortality rate.The mortality rate is a continuous function, Z(t), (11.2) and the density, f(t), is (11.3)If an age of first capture, t c , is considered, the simple model (4.17) leads to a battery of Weibull distributions corresponding to the unexploited phase subject only to natural mortality (with a coefficient of mortality M) and the exploited phase subject to total mortality (with a larger mortality coefficient, Z=M+F): If the stock process verifies condition (6.4), the corresponding stock equation (6.1) is (11.5)Under this model, when t ≥ t c , the exploitation rate (11.6) is constant.Therefore, the catch equation (8.2) yields (11.7)Thus, equations (11.5) and (11.7) generalize the classical equations of marine population dynamics which correspond to the deterministic version of the particular case of an exponent α =1.

THE SIMULATION OF SURVIVAL TIMES (DURATION OF LIFE)
Simulation procedures are a useful tool in both theoretical and applied research.They are dealt with in the general framework and in the particular proposed extended model.
Let T be a random survival time (life time) with any admissible survival, S(t), distribution function, F(t) = 1 -S(t), total mortality, Z(t), and cumulative hazard, H(t).
The classical simulation procedure is to obtain a sample of random numbers of a variable X which follows a uniform distribution in the unit interval [0,1] and to simulate the required lifetimes by applying the inverse of the distribution function: As an alternative, simple algebraic operations show that the random transformed lifetime H(t) follows a simple exponential distribution with a constant mortality rate equal to one.
Indeed, for the transformed lifetime, the survival is .
Using equation (1.5), S(t) = e -H(t) , and the decreasing character of the survival, we have where S -1 is the inverse of the survival function considered.
And e -t is the survival function of a standard exponential distribution of unit and constant mortality rate.
Let Exp(1) denote this exponential distribution, whose simulation is a standard tool in every statistical package, and let {x i , i = 1…n} be a simulated sample of size n of this distribution.The lemma implies that the set {H -1 , i = 1…n}, where H -1 is the inverse of the given cumulative hazard function, H, is a sequence of simulated lifetimes following the actual assumed model whose survival function is S(t).

S t T t S t S t e
As an example, the simulation of the extended proposed model given by the series of Weibull distributions (11.4) should be: (11.1) 13.REFERENCE PARAMETERS This section does not deal exhaustively with the biological and fishery reference point parameters in the specialised monographs (Caddy and Mahon, 1996;Cadima, 2000), but it analyses the most common reference points for management purposes.
Firstly, the principal tools generating the usual reference parameters were adapted to the proposed Weibull model.This is the case of biomass per recruit and yield per recruit, which generate some of the important reference points used for management purposes, such as the F MSY , F 0.1 , F med .They also provide important and useful concepts such as virgin biomass and overexploitation growth.For this adaptation, it was necessary to previously adapt the critical age.Secondly, we analysed some indices broadly used in all population dynamics (including human populations) but only marginally dealt with in fishery science, such as life expectancy, mean residual lifetime and median survival time.These parameters are redundant with mortality rates in the classical exponential model, but are not so trivial in a more general framework.The overall natural, fishing and total mortality rates were generalised.Finally, as a consequence of this, the relationship between an assumed level of natural mortality rate and the corresponding natural mortality coefficient in order to handle further extensions of the model considered in another paper (Ferrandis and Hernández, 2007).

The critical age
This reference parameter is defined as the age of maximum expected biomass in the evolution of the cohort.This biomass is (13.1.1)where: B(t) is the biomass of the cohort at age t. w i (t) is the weight of each individual at age t.Assuming the independence between the weight at any age and the survival where: N 0 is the initial size of the cohort, W ∞ is the asymptotic weight, K and t 0 are the usual Von Bertalanffy growth parameters, and b is the exponent of the weight-length relationship, then the critical age is the age that maximises the function whose derivative is Considering the definition of the total mortality rate, Z(t) = -S´(t)/S(t), this derivative may be expressed as (13.1.2) and equating to zero, it yields This implies that (13.1.3)which is the equation whose solution(s) give(s) the critical age(s) for any model Z(t) of the mortality rate.The practical computations are included in Appendix 1.

Life expectancy, mean residual lifetime and median survival time
Life expectancy and median survival time are parameters for describing survival time that can be useful in order to compare the different cohorts in a study period or the state of resources in different areas.
If T denotes the survival time whose survival function is S(t), the life expectancy, or the expectation of survival time (Smith, 2002), is given by: The usual statistical packages (SPSS, Statistika, S-Plus) have the means to calculate this cumulative distribution.
In the classical approach with a unit exponent, α = 1, the gamma function is one, the gamma distribution becomes the exponential one and the life expectancy is (13.2.3) that in the simple cases of an unexploited (t c = ∞) or a fully exploited resource (t c = 0) becomes, respectively, E(T) = 1/M and E(T) = 1/Z.In the author's opinion this trivial redundancy of the life expectancy with the mortality rates in the classical approach is the reason why the life expectancy, considered as a fundamental reference parameter in general population dynamics (including that of human populations), is systematically neglected in fishery science.
An interesting parameter related to the life expectancy is the mean residual lifetime at age t, r(t), as a measure of the average remaining life.This is given by (13.2.4) where the integral is calculated in the same way as is applied for the life expectancy (Dobson, 2002).) .( / , ) In particular, for the age of first capture, t c .
Notice that in the particular and classical case of a unit exponent, α = 1, the mean residual time of life should be r(t c ) = 1/Z.
The median survival time, t 50 , is given by the solution to the equation and is often considered (Dobson, 2002) as a better description of the "average" survival time than the life expectancy because of the bias of the survival distribution.
The median is a particular case of a percentile, t p , of level p of the lifetime distribution, i.e. when p = 50%.Thus In the current case of a series of Weibull distributions given by (11.4), we have 13.3.The mortality rates

The natural mortality rate
The constant natural mortality rates established in the literature and estimated by different methods can be adapted to any admissible mortality model, and particularly to the proposed Weibull model.A constant mortality rate M 0 corresponds to the assumption of an unexploited population with an exponentially distributed lifetime and a life expectancy E(t) = 1/M 0 .
For the Weibull composite model, setting t c = ∞ in (13.2.2), we have (13.3.1) and (13.3.2) Considering, as is assumed in the proposed survival model (11.4), the exponent of the Weibull distribution as a population parameter and the life expectancy as an overall parameter of the lifetime distribution, the Weibull distribution with parameters M 1 and α would have the same life expectancy as that attributed to the virgin population, affected only by natural mortality.
The overall or expected mortality rate, E(M(t)), corresponding to a given natural mortality coefficient, M, for an unexploited population is which, with the change of variable, x = Mt α , leads to (13.3.3)Now, taking the constant assumed natural mortality rate M 0 as an estimation of this expected (or mean) mortality rate, i.e.M 0 = M -[0, ∞), we have an alternative expression for the coefficient of natural mortality: (13.3.4) The Weibull distribution with parameters M 2 and α would have the same overall expected mortality rate as the one attributed to the virgin population, which is affected only by the assumed and constant natural mortality.
Comparing the two estimates of the coefficient of natural mortality, their ratio is which, for α > 0.5 is ≥ 1, by the properties of the gamma function.Thus, the estimate obtained by standardising the life expectancy (Expression (13.3.2)) will be generally greater.

M M
For the particular case of the usual exponential distribution, we have Thus, the two estimates are the same and equal to the assumed M 0 .

The overall fishing mortality rate
In the proposed model, the mortality rates vary with age.The overall (expected or mean) fishing mortality rate in the interval [t c , ∞), from the effective age of first capture, t c, is given by .
For the proposed Weibull model, , which again, with the change of variable x = Zt α , yields , which gives (13.3.5)

The overall total mortality rate
The overall total mortality rate is given by .Again, for the proposed Weibull model, A similar development as the above one gives (13.3.6)

The overall exploitation rate
The overall exploitation rate is given by (13.3.7) In the case of the proposed Weibull model, we have

The yield per recruit
The yield per recruit, Y, is a reference parameter frequently used for diagnosis of the actual effort in relation to the maximum sustainable yield, MSY, and its precautionary limits.It is defined as the relative (by recruit) biomass catch from a given cohort. Its or, in terms of the conditional survival from the age t c (13.4.2) In order to calculate the integral with a given precision, ε, the infinite interval [t c ,∞], is split into two parts: a finite interval [t c ,n] bounded by an integer n, where the approximation to the integral may be obtained, and an infinite interval [n,∞), where the value of the integral should be negligible.The concrete development is presented in Appendix 2.

The biomass per recruit
As in the previous paragraph, t r let designate the age of recruitment, and consider the age interval since recruitment, [t r ,∞), split into the two intervals [t r ,t c ) and [t c ,∞).
Its differential equation is Hence, the expected biomass per recruit is given by (13.5.1) If the fishing mortality increases, the survival function decreases and the biomass per recruit is obviously decreasing with respect to the fishing mortality, as in the classical approach by Beverton and Holt.Both integrals are particular cases of the computational tools described in Appendix 2 for the case of the yield per recruit.

CONCLUSIONS
The relation between survival, mortality and lifetime distribution is the basis on which a stochastic approach to marine population dynamics is built.The properties of survival and mortality have been rigorously established and examples of published and unpublished, adequate (admissible) and inadequate models have been made.Specifically, alternative and generalised models using the Beverton and Holt and Sparholt hypotheses on natural early mortality have been presented.They are admissible mortality models continuously decreasing with age according with the ideas elaborated by Caddy (1991).
The analysis of total mortality rate in the two independent components, natural and fishing, is essential in the present elaboration as it is in all marine population dynamics.However, in the future it may give rise to dependent causes of death or "competitive risks", a concept that may be useful in an ecosystem approach framework.
The definition of survival provides an initial binomial stochastic process for the evolution of the stock of a cohort of marine resources and the population growth model leads to a lognormal and more complex version of the process.Both are Markovian processes providing the fundamental equation of population dynamics as a conditional expectation.In the present work, both binomial and lognormal processes are treated in a unified approach.These processes are correctly defined with their Gaussian approximations, which generate a promising bridge between marine population dynamics and multivariate statistical analysis.The identification of the determinant and the inverse of the covariance matrix converts this Gaussian processes into a "ready to use" tool.
The catch process is also correctly defined.It provides the catch equation as a conditional expectation that is easily formulated when the mortalities corresponding to natural, fishing and total mortality are proportional.This condition is related to the Weibull distribution for the lifetime of the population.
Thus, an extension of the classical marine population dynamics is proposed through survival and lifetime models as a battery of two Weibull distributions, considering the age of first capture as the cutoff of the unexploited and exploited phases of the cohort.
Under this extended model, the reference parameters related to the diagnosis of the fishery have been established in a way that generalises the known expressions derived from the classical exponential model.
The critical age, life expectancy, median survival time, mean residual lifetime, natural, fishing and total mortality rates, biomass per recruit and yield per recruit have been generalised and computer solutions have been implemented.
The idea underlying the Weibull survival model here presented is an alternative to classical population dynamics.Instead of considering many age intervals with constant mortalities, i.e. the mortalities as step functions, we can establish few intervals but with more flexible mortality models, i.e. the mortalities as continuous functions.
According with the author's experience with Mediterranean demersal resources, this Weibull survival model has a strong coherence with the estimated survival functions of some target species.It is for this reason that the hypotheses assumed by this model are not significantly contradictory with the behaviour of trawl surveys and commercial data.This does not hold for a general application of the proposed Weibull model.In a few words, the Weibull model may be reasonable and adequate for some fishery resources, though neither perfect nor generally applicable.The elaboration of the Beverton and Holt and Sparholt mortality models here presented shows a way (not the only one, nor the best) to deal with more complex situations.

APPENDIX 1. THE CRITICAL AGE
In the case of a constant mortality rate Z(t) = Z, Equation (13.1.3)gives the unique and well-known solution for the critical age.
In the more flexible case of an exponent α ≠ 1 corresponding to the Weibull model proposed, by the properties of the survival function formalised in the first paragraph and thus, for any exponent α < 0 because .
First, the alternative α > 1 will be considered.Now, the mortality rate, Z(t), is a monotone increasing function Then, the second factor of Expression (13.1.2) is a monotone decreasing function which approaches -∞ and which has a unique discontinuity at the age of first capture t c .
Therefore, the derivative of the expected biomass of the cohort (Expression (13.1.2))has the following properties: , because t 0 < 0; Let us now consider the case when the exponent is lower than one, 0 < α < 1, leading to a piecewise continuous and decreasing mortality rate.
We can rewrite Expression (13.1.2) in order to identify its asymptotic sign which will depend on the product Z(t)e  Consequently, if B´(t) reaches any positive value, the range of ages [0,∞) will contain an interval [t cr ,∞), in which B´(t) will have only negative values that change to positive below the cut-off age, t cr .This age will therefore correspond to a relative maximum of the biomass of the cohort, and hence to the critical age.
As lim t→∞ B´(t) = -∞ , an initial interval [0, t n ) will exist in which B´(t) will also present negative values that change to positive beyond the cut-off age, t n , which will correspond to a relative minimum of the biomass of the cohort.
A sufficient and easy condition for verifying the existence of positive values of the derivative, B´(t), is to check whether B´(t c ) > 0. This inequality implies condition (A.1.2)and, in this case, the critical age will be in the exploited phase, t > t c , and the solution to (A.1.3).
Solutions to (A.1.1)or (A.1.3)must be obtained by computer numerical methods like the false position (regula falsi) method (Hildebrand, 1974;Curtis and Patrick, 1994).The author has established suitable software for solving these equations with a prefixed precision.APPENDIX 2. THE YIELD PER RECRUIT where The weight curve reaches its inflexion point at the age .Hence, the second term is and due to the integrability of the survival , for any given ε > 0 it is possible to choose n such that B ≥ ε/2.In fact, for the proposed composite Weibull model but and then, it is possible to choose an integer n such that Once an n is identified, the approximation of the first term, A, may be easily obtained by the trapezoidal or Simpson methods.
The author has built specific programs for calculating the yield per recruit with a given precision, ε, prefixed by the user.
corresponds to the original Beverton and Holt relation between mortality and survivors: (4.10) conditional survival for the individuals who reach the age t 1 is function.The integral is the cumulative distribution function of the standard gamma distribution of shape parameter p. Then and, finally, the life expectancy is (13.2.2) ., 71(1),March 2007, 153-174.ISSN: 0214-8358 ., 71(1), March 2007, 153-174.ISSN: 0214-8358 STOCHASTIC POPULATION DYNAMICS • 173 because Therefore, B´(t), from a certain age it will present only negative values.Now we can summarise the properties of the derivative of the expected biomass of the cohort: B´(t) has a unique discontinuity at the age of first capture t c ; B´(t) reaches and retains negative values from a certain age.
has a unique discontinuity at the age of first capture t c ; B´(t) reaches and retains negative values from a certain age (the critical age) .Hence, in this case, a unique critical age exists.It will be in the unexploited phase, t < t c , and the solution of the critical age will be the age of first capture t c .