Aqua-BioSci. Monogr. (ABSM), Vol. 2, No. 3, pp. 1-45 (2009)

doi:10.5047/absm.2009.00203.0001

www.terrapub.co.jp/onlinemonographs/absm/

# Non-linear and Graphical Methods for Fish Stock Analysis with Statistical Modeling

## Tatsuro Akamine*

National Research Institute of Fisheries Science, Fisheries Research Agency,

Fukuura 2-12-4, Yokohama 236-8648, Japan

### Abstract

Useful methods for growth curve fitting, body-size composition analysis, and estimation of population size in fish stocks are presented. These methods are statistically based on the maximum likelihood method and the likelihood ratio test. Mathematical explanation of the standard Richards growth formula with seasonal change, the generalized reproduction model, and the Awaya method for estimating implicit function models are given. Mathematical proofs of the iteration method, called the Hasselblad method, or the EM algorithm for estimating the mixture of normal distributions, and the Marquardt method for general optimization are shown. For population size estimation, the Petersen method for mark–recapture experiments, the quadrat method, and the DeLury removal method are discussed. These are based on the binomial distribution and the classical Bayesian statistical methods which are also discussed. Mathematical proofs of the sum formulae of the binomial and hyper-geometric distributions are given. The virtual population analysis using mortality rates, the Leslie matrix model, and the linear programming for discrete fishing models are also explained. All the methods stated here can be easily carried out using spread-sheet software.

### Keywords

body-size composition, growth curve, population size, reproduction, survival

Received on March 4, 2008

Revised on March 13, 2009

Accepted on April 28, 2009

Published online on July 28, 2009

*Corresponding author at:

National Research Institute of Fisheries Science, Fisheries Research Agency, Fukuura 2-12-4, Yokohama 236-8648, Japan

e-mail: akabe@affrc.go.jp

### Introduction

There are many methods to estimate the parameters of fish stocks. Before the 1980s, most of them were approximations, linear regression methods, and generally had low precision. Since the 1980s, the author has developed various non-linear methods by using microcomputers. Nowadays, these methods can be carried out using spread-sheet software. In particular, its graphical functions are useful to understand the statistical models and check parameter values.

In Section 2, the standard growth formula of fish is presented. This formula is based on the Richards growth formula and expanded for seasonal growth change. The non-linear regression method to estimate and test parameters, the generalized reproduction model, and the Awaya method for estimating implicit function models are also explained. In Section 3, the Hasselblad method which estimates parameters of a mixture of normal distributions is explained by using the EM algorithm and the K–L information quantity. The proof of convergence, the Marquardt method and the Jacobi method are also shown. In Section 4, the estimation methods for population size are detailed. These are the Petersen method for mark–recapture experiments, the quadrat method, and the DeLury removal method. Their statistical models are based on the binomial or hyper-geometric distribution. The classical Bayesian statistical methods and the sum formulae of the binomial or hyper-geometric distributions are explained. Survival models are shown in Section 5. The virtual population analysis (VPA) using mortality rates, the Leslie matrix model related to VPA, and the linear programming method for discrete fishing models are explained.

The above-mentioned models are simple and easy to understand. Statistical methods of point and interval estimations are the maximum likelihood method and the likelihood ratio test. We can easily calculate and draw figures for statistical estimation by using spread-sheet software.

page top

### 2. Standard growth formula in fish population dynamics

In fish population dynamics, many growth formulae and their expansions are used. Among them, von Bertalanffy growth formula is the most commonly used. For parameter estimation, the linear regression method (e.g. Ford/Walford plot) for data measured with a constant interval, has previously been widely used. However, it was impossible to estimate three parameters at the same time and the curve fitting was not good. In the 1980s, personal computers and BASIC made it possible to estimate non-linear models and obtain good curve-fitting. In this section, some expanded models, the standard growth formula, and estimation method will be shown.

page top

### 2-1. Traditional growth formulae

In fish population dynamics, the von Bertalanffy growth formulae are widely used for body length (VBGF1)

and for body weight (VBGF2)

where l and w are the body-sizes, t is the time, k is the growth coefficient, l and w are the body-sizes at t = ∞, t0 is the extrapolated time when l = 0 or w = 0. These are solutions of the differential equation

where η, κ and a are constants. On the other hand, the logistic growth formula (LGF)

and the Gompertz growth formula (GGF)

where c is the time of the inflection point, are also used in some cases.

Richards (1959) solved the differential equation

and obtained the growth formula (RGF)

or

with the initial condition (t, w) = (0, w0). This formula corresponds to LGF when m = 2, VBGF1 when m = 0, VBGF2 when m =2/3, and converges to GGF when m → 1, m ≠1. However, this convergence is not easy to understand because Eq. (2.6) is an exponential function when m = 1 (Taylor 1962).

In ecology and agriculture, France and Thornley (1984) defined RGF as follows:

where wf = w. This is the solution of the differential equation

However, parameters are different from the traditional growth formulae used in fish population dynamics.

For use in fish population dynamics, Pauly and David (1981) suggested the following formula:

This is the generalized VBGF. Although Schnute (1981) defined RGF as

this is the generalized LGF which is incorporated in RGF. He also defined the generalized VBGF

and GGF

These four equations are included in RGF. The application and potential of RGF has not been sufficiently understood in fish population dynamics.

page top

### 2-2. Standard formula of RGF in fish population dynamics

Akamine (1988b) defined the differential equation of RGF as follows:

where r is the parameter which decides the shape of the growth curve. It is easy to understand the convergence when r → 0 as follows:

by using the definition of the logarithm

Equation (2.16) is the differential equation of GGF. The transformation

reduces Eq. (2.15) to the following differential equation

The solution of this with the initial condition (t, w) = (c, w/(1 + r)1/r) is

This is the standard formula of RGF used in fish population dynamics, which also includes generalized VBGF, GGF and generalized LGF (Fig. 1). The initial condition means the inflection point of this curve. We can understand that Eq. (2.20) converges to Eq. (2.5) when r → 0 by using the following definition of the exponential:

When r = −p < 0, this standard formula can be rewritten as

This is the generalized VBGF. By using the allometric function

RGF is transformed to

This is also RGF.

Fig. 1. Graphs of the Richards growth formulae. Figures are values of r. Revised from Bull. Japan Sea Natl. Fish. Res. Inst., 38, Fig. 2, Akamine T. Estimation of parameters for Richards model. 187–200, 1988, with permission from Japan Sea National Fisheries Research Institute, Fisheries Research Agency.

Schnute (1981) defined his new growth formula as the differential equations

Equation (2.15) can be rewritten as

These are modified as follows:

Therefore, Eq. (2.25) is equal to Eq. (2.15) which is RGF itself.

page top

### 2-3. Seasonal growth formula

Pitcher and MacDonald (1973) presented the growth formulae (PM1)

and (PM2)

where tr is real time (weeks), tg is growth time, ts = trs, s is the starting point for cosine, sw is the threshold value for cosine, C is the magnitude constant, and s1 is the starting point for sine. PM1 is the switching model and PM2 is the sine curve model. And they define the water temperature θ as

where m is the mean water temperature, q is the magnitude constatnt, and s2 is the time lag between temperature and growth switch. Using this function, they rewrote the growth formulae as follows (PM1):

and (PM2):

where sw1 is the threshold temperature, and s3 is the time lag between mean temperature and growth oscilation. Pauly and Gashutz (1979) presented the equation (PG)

This is simpler than PM2. Haddon (2001) used the formula (HD)

This is a more complex model than PM2 and PG. Appeldoorn (1987) suggested that these formulae have a problem as L(t0) ≠ 0.

Cloern and Nichols (1978) defined the differential equation

where L'max is maximal growth rate when body size is minimal, and Lmax is the upper limit of body size. They obtained the solution (CN)

Although this model overcomes the problem stated by Appeldoorn, this solution has two further problems. The first is that the parameter Lmin is not necessary and the second is that the amplitude of the periodic function is fixed to be 1.

Hoenig and Hanumara presented the formula in 1982 (HH)

Furthermore, Somers (1988) presented the formula (SO)

This is a simpler version of HH. These models also overcome the problem stated by Appeldoorn. In summary, researchers in fish population dynamics have previously used the following expression for exponential:

However, the meaning of the trigonometric function f(t) is difficult to understand.

page top

### 2-4. Standard formula for seasonal growth

Akamine (1986) solved the differential equation

and obtained the growth formula (AK1)

This is the simplest expression for seasonal growth of VBGF1. He used the functions

Where af(t) ≤ 1. In the case that f(t) is the water temperature, F(t) means the cumulative water temperature. When a < 0, this formula includes negative growth. He also obtained an expanded expression of the logistic growth formula (AK2)

and that of the Gompertz growth formula (AK3)

Further, Akamine (1988b) obtained an expanded Richards growth formula (RA):

This is the standard formula for seasonal growth. Akamine (1993, 1995) showed the general case as follows: Let the differential equation be

He later modified this to

The integral of this is

where C is an integral constant. Then he obtained the general solution

and the particular solution

Therefore, we can get the expanded formula with the operation tF(t).

He used the following standardization:

and functions

Equation (2.43) is rewritten as

In this formula, the parameter a strongly influences the growth coefficient k. Therefore, the following operations are needed for comparison of parameter values:

Kiso et al. (1992) applied the following function to data for the Masu salmon (Oncorhynchus masou):

The standard formula (2.46) can include many growth formulae. For example, it includes the switching model PM1. Let

then we obtain

This is equal to PM1. Furthermore, Pauly et al. (1992) proposed a new switching model (PE):

where NGT is the duration of non-growth within a year. This is also included in the standard formula (2.46). Let

then we obtain

This is equal to PE.

page top

### 2-5. Parameter estimation and statistical test of growth formulae

For parameter estimation of growth formulae, Pitcher and MacDonald (1973) calculated the fitting by hand and the direct-search method by computer. They also used the multivariable regression to obtain values of C, k and t0 in Eq. (2.29) as follows:

Pauly and Gashutz (1979) also used the multiple regression to obtain values of C, k, t0 and ts in Eq. (2.33) as follows:

These are linear regression models that cannot obtain a high precision estimate for L.

Schnute (1981) used the simplex search method, which was transformed to microcomputers, for his growth model. The objective functions are the least-squares method with no weighting for raw data, and for their logarithm values.

Akamine (1982) used the Gauss–Seidel method in BASIC to estimate parameters of a mixture of normal distributions and Akamine (1984) used the Marquardt method for the same model. These objective functions are the least-squares method with no weighting. Akamine (1986) used the Marquardt method for the estimation of the growth curves and the objective function is the weighted least-squares method.

When we measure fish body-size w at time ti, let ni be the number of specimens, be the mean, and σ2i be the variance. When t = ti the residual sum of squares is

where the subscript i is omitted, wj is the body-size of the j-th individual at time ti, and w(t) is the growth formula. Therefore, we can calculate the residual sum of squares only by using the mean and its variance. According to the central limit theorem, the distribution of is the normal distribution and its variance is given

Thus, the objective function of the weighted least-squares method is defined as

The optimization method which minimizes this objective function gives the parameter values of the growth formula.

(Example 1) Fitting the growth formula to clam data.

The clam data of the Oita prefecture is shown in Table 1. Although these data have the means and their standard deviations, they have no number of specimens. Therefore, let δ ≡ σ and use the weighted least-squares method (2.67). By using the optimization method of the spread-sheet software (Gorie 2001), we obtain the Gompertz formula

From this solution as the initial values, we obtain the expanded growth formula for seasonal growth

Because A > 1, this growth formula has periods of negative growth. This is not valid for hard structures such as shells. These data are for samples taken only four times in a year. It requires more samples for a better estimation. Figure 2 is given by the graphical function of spread-sheet software. This curve has periods of little negative growth because software interpolates. It seems more natural for the real growth pattern than the exact drawing.

Table 1. The growth data of the clam (Meretrix lusoria). SL is the Shell length (mm). Data were excerpted from Kamijoh et al. Hamaguri no shigen baiyou gijutu kaihatu kenkyu houkokusho, 1985, with permission from Fisheries Research Institute, Oita Prefectural Agriculture, Forestry and Fisheries Research Center.

Fig. 2. The growth curve of the clam. Black circles are SL, and squares are S.D. Revised with permission from Akamine T. Taicho-sosei data karano nenrei to seicho no suitei. In: Shigen Hyoka notameno Suuchi Kaiseki, Suisangaku Series 66, p. 127, Fig. 8.6, Kouseisha Kouseikaku, Tokyo, 1987. © 1987, the Japanese Society of Fisheries Science.

The interval estimation and the test of the growth formula are easy because the distribution of Y is the χ2 distribution (Akamine 2004). The null hypothesis

where p is the number of parameters to estimate, gives the following distribution:

where Y0 is the objective function of the true values and Ymin is that of the point estimate. We can obtain the confidence area of parameters by using this distribution. On the other hand, the AIC (Akaike Information Criterion) is defined as

where ln Lmax is the maximum log likelihood. In this case, the AIC is rewritten as

When there is no weight, the minimum of the residual squares

where m is the number of specimens, gives the optimum parameter values. As same as for the linear regression, we can use the F distribution for the interval estimation or the test as follows:

where S0 is the residual squares of the true values and Smin is that of the point estimate. On the other hand, the AIC is approximated in this case as

The example of the model selection by using the AIC is shown in Kiso et al. (1992).

As a typical case, the test of the growth difference between male and female is shown as follows: Let the number of male data be mM, that of female data be mF. The null hypothesis

is used for this test. For the weighted case, we can use the following distribution:

On the other hand, for the non-weighted case, we can use

Although the confidence area of the linear model is an ellipse, that of the non-linear one is a banana-shape. In fisheries science, the wrong discussion has sometimes been seen in which the overlapping between male and female confidence intervals means that the same growth curve can be fitted. Because the confidence interval does not mean the probability in which the true value exists, it is better to use the null hypothesis for these subjects.

(Example 2) Fitting VBGF1 for Pacific hake data.

The data is shown in Table 2. This data has no variance. Thus, F test (2.77) is adequate. The optimization method of the spread-sheet software gives two VBGF1s for male and female by using Eq. (2.72) as follows (Fig. 3):

On the other hand, the total data of male and female gives

Thus, Eq. (2.77) is

This is significant at 5%, but not at 1%, because F(3,18) = 3.160 (5%), 5.902 (1%).

Table 2. The body-length data (cm) of the Pacific hake (Merluccius productus). Data were excerpted with permission from Fishery Bulletin, 73, Dark TA. Age and growth of Pacific hake, Merluccius productus. 336–355, 1975.

Fig. 3. Black circles and the thick line are female data and white circles and the thin line are male. Redrawn with permission after Fishery Bulletin, 77, Fig. 1, Kimura DK. Likelihood methods for the von Bertalanffy growth curve. 765–775, 1980.

The confidence region of female growth parameters can be estimated using Eq. (2.73), which is rewritten as

We can obtain the confidence region of parameters by using S0. The null hypothesis must be

In this case p = 3, m = 13, F(3,10) = 3.708 (5%), Smin = SF = 28.8, and we obtain S0 = 60.9 for the condition. By using the optimization method of spread-sheet software, we obtain the minimum k = 0.200 (l = 66.3, t0 = −0.761), and the maximum k = 0.413 (l = 58.3, t0 = 0.346), where these values are edges of the three-dimensional confidence region. These VBGF1s are shown as thick lines in Fig. 4. In general, there is a positive correlation between k and t0, negative correlations between k and l, and l and t0 in VBGF1. The minimum and maximum estimates of l and t0 are omitted.

Fig. 4. The 95% confidence interval of the growth curve for female Pacific hake. Two thick lines show growth curves for the minimum and maximum k. Redrawn with permission after Akamine T, Suisan Shigen Kaiseki no Kiso, p. 32, Fig. 2.6. Kouseisha Kouseikaku, Tokyo, 2007. © 2007, Kouseisha Kouseikaku Co., Ltd.

page top

### 2-6. The generalized reproduction model

The treatment of the reproduction model is similar to that of the growth formula. A generalized reproduction model was proposed by Schnute (1985):

This model can be rewritten as

where S is the spawning biomass, R is the recruit biomass, and α, β and r are the parameters which determine the shape of the curve. When r = 1, Eq. (2.81) is equal to the Beverton–Holt model

When r → 0, Eq. (2.81) converges to the Ricker model

When r = −1, Eq. (2.81) is equal to the parabola which is the Shaeffer model of the surplus production model. When r → ∞, Eq. (2.81) converges to a straight line (Fig. 5). The gradient at the origin is α. These curves have the following characteristics:

1. When r < 0, they have an intersection point or a contact point on the abscissa at S= −1/(rβ). The right side over an intersection point or a contact point of these reproduction models has no meaning.
2. When r < 1, they have a maximal point at S = 1/[(1 − r)β].
3. When −1 < r < 1, they have an inflection point at S = 2/[(1 − r)β].

Fig. 5. The graphs of the generalized reproduction model. Figures are values of r. Redrawn with permission after Schnute J. A general theory for analysis of catch and effort data. Can. J. Fish. Aquat. Sci., 42(3), 419–429, Fig. 2, 1985. © 2008, NRC Canada.

Formerly, we obtain this model as follows: The survival equation is

where N is the number of fish. The solution of this equation with the initial condition N(0) = S is

Substitution of R = N(t), α = et and β = b(1 − ert)/r into this equation gives Eq. (2.81). On the other hand, the following differential equation gives the same solution (Akamine 1994):

The transformation

gives the equation

The solution of this with the initial condition R'(0) = α is Eq. (2.81).

page top

### 2-7. Parameter estimation for reproduction model

In the case of estimating R from S, it is better to generally use the regression theory. If we consider both errors of R and S, the weighted least-squares method with errors in both axes is applicable as follows: Let the reproduction model be the implicit model

The objective function with m-paired data (si, ri, δsi2, δri2):

is adequate for this estimation (Awaya 1991). We can minimize Eq. (2.90) by using Lagrange's indetermined multiplier method.

Awaya (1983) applied the Newton method for the implicit function model. Let the implicit function be

where t is an independent variable, w is a subordinate variable, and θ are parameters. This function cannot be rewritten as w = g(t, θ). Let the data be m sets of (ti, μii2) and the objective function be

We can obtain adequate parameter values by minimizing this objective function. Let the objective function be expanded as follows:

where vi = 1/δi2, and λi is the undetermined multiplier. Thus, the following simultaneous equations can be obtained:

where fi = f(ti, wi, θ), and n is the number of parameters. Let

where the subscript 0 means the initial values. The linear approximation of Eq. (2.96) is

where fi0 = f(ti, wi0, θ0). Thus, the following equations can be obtained:

On the other hand, substituting Eq. (2.97) into Eq. (2.94) with fi = fi0 gives

where di0 = μiwi0. Substituting Eq. (2.100) into Eq. (2.101) gives

Substituting these into Eq. (2.95) with fi = fi0 gives the following equations:

These are the normal equations of the least-squares method for the implicit function model.

In practice, if we first decide the initial values of θi, the initial values of wi, which satisfy Eq. (2.96), can be obtained by using numerical calculation. Solving Eq. (2.103) for the modified vector of parameters Δθi gives the modified vector of variables Δwi from Eq. (2.100) and new values of fi from Eq. (2.99). Therefore, solving Eq. (2.103) again and repeating these procedures gives the solution. The linear approximation errors of fi converge to 0. This Awaya method is regarded as the Newton method for the multivariable implicit function model. Although its convergence area is narrow, it converges in 10 times or so. This method can be modified to the Marquardt method which has a larger convergence area (Awaya 1991; Akamine 1999). The merit of the Awaya method is that it requires only one time to calculate the correct value of wi.

page top

### 2-8. Supplement

In fish population dynamics, allometry

is used for the relationship between the body-weight and body-length. Nowadays, we can estimate body-weight from body-length data by using the following two regression models:

where ε is the residual. Each objective function is

Y2 is better than Y1 because Y2 is expected to have the same variance for each data. On the other hand, the objective function (2.90) is better to estimate parameter values a and b when both w and l have errors. The objective function must be selected for the purpose.

Ohnishi and Akamine (2006) presented the growth formula for shells (OA). The differential equation of this is

where s = ξ l2, w = ηl3. Thus we obtain the equation

OA is equal to VBGF1 when a = 0 and converges to LGF when a → ∞. Although OA is similar to RGF, it does not include GGF and VBGF2. The integral of OA is

where A = a/l, l0 = l(0). This is an implicit model, so the Awaya method is useful for parameter estimation.

page top

### 3. Analysis of the body-size composition

Estimating the age composition from the body-size composition data is a basic method in fish population analysis. Before the 1960s, graphical methods for a mixture of normal distributions were widely used. Hasselblad (1966) invented the iteration method for this model by using the maximum likelihood method. Akamine (1982, 1984, 1985, 1987b, 1988a) proposed BASIC programs for the Gauss–Seidel method and the Marquardt method by using the maximum likelihood method and the least-squares method. By use of spread-sheet software, the Hasselblad method is easy to apply. In this section, the Hasselblad method is explained by using the EM algorithm. The Jacobi method and the Marquardt method are also explained.

page top

### 3-1. Statistical model

The body-size composition is regarded as a mixture of normal distributions which is defined as

where x is the body-size, subscript i is the number of each distribution, and pi > 0 is a mixture ratio which has the condition

Function gi(x) ≥ 0 is a normal distribution:

The estimates of pi, μi and σi2 can be statistically obtained from the body-size composition data. The condition (3.2) decreases the freedom of parameters. Thus, we have to estimate (3n − 1) parameters at the same time.

The maximum likelihood method is adequate to estimate parameters of the probability distribution. When the body-size data is obtained as x1, ..., xm, the likelihood for this model is defined as

The parameter values which give the maximum likelihood are the most adequate estimation. In practice, the log likelihood is used for the objective function as follows:

In the case of the histogram for the data, let H*(x) be the frequency of a histogram, and H(x) > 0 be the relative frequency

where x is the class mark of the body size. Thus, the log likelihood is defined as

We can obtain the maximum Y by using the iteration method for computers.

page top

### 3-2. Hasselblad method

Hasselblad (1966) defined the mixture ratio as follows:

In this case, the condition of the solution

is rewritten as

where

By multiplying pi to both sides of Eq. (3.10), he obtained

The sum of both sides are

Therefore,

The resulting equation from Eqs. (3.12) and (3.15) is

Thus, he obtained the following iteration method:

The right upper (t) means the t-th iteration value. Substituting old values in the right side provides new values in the left side.

On the other hand, the condition of the solution:

is rewritten as

Thus, he obtained the following iteration:

In the same way, the condition of the solution

is rewritten as

Thus, he obtained the iteration

Although the Hasselblad method is stable and useful, it is difficult to fully understand the properties of iteration (3.17).

page top

### 3-3. Undetermined multiplier method

Akamine and Matsumiya (1992) obtained Eq. (3.17) by using Lagrange's undetermined multiplier method. Let the objective function be

where λ is an undetermined multiplier. The condition of the solution

is rewritten as

The value of λ is determined as

Thus, we obtain the following non-linear equations:

By multiplying pi to both sides of this equation, we obtain Eq. (3.16) and iteration (3.17).

page top

### 3-4. EM algorithm

The Hasselblad method is regarded as an application of the EM algorithm which has E-step for estimating the missing data and M-step for maximizing the likelihood to estimate parameters. The EM algorithm is developed for the model which has random missing data. Although the model of a mixture of normal distributions has no missing data, it has the latent data as follows:

Define the latent data zi(x) as the relative frequency of age i and size x, and qi(x) as the probability of age i in the size x as follows:

In E-step, zi(x) can be estimated as the expectation:

These equations are natural and the calculation amount of E-step is zero.

For the optimization of parameters in M-step, there are defined relations

We obtain the following log likelihood:

To maximize this Q-function for the parameter pi, let the objective function be

The condition of the solution

is rewritten as

The value of λ is calculated as

Thus, we obtain the iteration

This is the same as Eq. (3.17). For the parameter μi and σi2, we also obtain the same iteration as the Hasselblad method. Therefore, we can obtain the iterations naturally by maximizing the Q-function in the EM algorithm.

It is proved that the log likelihood increases in each step of the EM algorithm. Akamine (2007) used the geometrical method (Murata and Ikeda 2000) for this proof. For discrete models, the K–L information quantity is defined as

where ξ(x) and η(x) are probability distributions. This inequality is proved by using ln xx − 1 as follows (Sakamoto et al. 1986):

The equality is clearly satisfied only when ξ(x) ≡ η(x). The K–L information quantity gives the upper limit of the log likelihood as

Define the distance as

This is a special distance because

Firstly, to minimize D we can modify zi(x) in E-step. We can divide D as follows:

Therefore, the equations

make D minimum because the first term of Eq. (3.43) is 0 and the second term does not include zi(x).

Secondly, to minimize D we can modify the parameters in M-step. The distance D is divided in another way

Because the first term is constant, we have to minimize the second term which is the same as the minus Q-function (3.32). Thus, it was proved that the distance D is made shorter in each E- and M-step. In E-step, D is modified as

In M-step, D is minimized. Thus, the objective function Y is increased in each step. It has an upper limit as Eq. (3.40). Then it must be converged. Although the convergence of the objective function Y is proved, the convergence of each parameter has not yet been proved.

page top

### 3-5. Convergence criterion by diminishing mapping

On the estimation of the mixture ratio, Akamine (2001) showed the following convergence criterion according to Iri (1981). Iteration (3.17) can be rewritten as

We obtain the following inequality:

where |J(t)| is the norm of the matrix J(t). The elements of this are

where the ∑ with ik means the sum of i = 1, ..., n except i = k. Thus, we obtain the following criterion:

where L is the Lipschitz constant. This is the sufficient condition of the convergence. On the other hand, the necessary condition is

page top

### 3-6. Approximation of the Jacobi method

Although the Hasselblad method is regarded as the EM algorithm, as detailed above, it is also regarded as an approximation of the Jacobi method in which all parameters are modified by using the Newton method at the same time in each iteration (Akamine 1998, 1999). Let Eq. (3.28) be

Thus, we obtain the following equation:

In the neighborhood of the solution, we obtain the approximation

where Ψi is Eq. (3.50). Therefore, the modification of pi in the Jacobi method is

This shows that the Hasselblad method is regarded as an approximation of the Jacobi method when Ψi is small.

page top

### 3-7. Difference between the iteration method and the EM algorithm

There is little difference between the iteration method and the EM algorithm (Akamine 2005). It is clear when σi = αμi. In this case of the EM algorithm, the condition of the solution

gives the equation

Therefore, we obtain the iteration

On the other hand, the condition of the solution

gives the quadratic equation

The iteration is

Equations (3.59) and (3.62) are impossible to solve because the values of α and μ on the right side are new values. In the iteration method, we can solve these equations by making the values of α or μ on the right side be old values, which is equal to the Hasselblad method. Therefore, in this case, M-step of the EM algorithm needs the iteration method. On the other hand, Akamine (2005) presented the Newton method for Eq. (3.61) as

This method converges in the same way as the Hasselblad method.

(Example 3) Estimation of the age composition for the data of the porgy.

The body-length composition is shown in Table 3 and Fig. 6. Let the histogram data be H*(7.5) = 7, H*(8.5) = 79, and so on. The iterations (3.17), (3.20) and (3.23) can be calculated by using the spread-sheet software (Aizawa and Takiguchi 1999; Gorie 2002). Table 4 shows the results.

Table 3. The body-length composition data of the porgy (Dentex tumifrons). BL is the body length (cm). Data were excerpted with kind permission of Dr. Shoichi Tanaka from Table 2 in Tanaka (1956) with its original data source of Group of Statistics, Pelagic Resource Section, SEIKAI Regional Fisheries Research Laboratory, Data of the investigation of September to November in 1950, Report of the Investigations of Trawl Fishery in the China Sea, No. 2, pp. 1–28, 1951. SEIKAI Regional Fisheries Research Laboratory, Nagasaki, Japan.

Fig. 6. Result of the estimation (Tabe 4). Black circles are data, thick line is a mixture, and thin lines are normal distributions. Redrawn with permission after Akamine T. Suisan Shigen Kaiseki no Kiso, p. 6, Fig. 1.3, Kouseisha Kouseikaku, Tokyo, 2007. © 2007, Kouseisha Kouseikaku Co., Ltd.

Table 4. The result of the estimation for the mixture of normal distributions. The range of the body length is 7–32 cm. Data were excerpted with permission from Akamine T. Suisan Shigen Kaiseki no Kiso, Table 1.2, Kouseisha Kouseikaku, Tokyo, 2007. © 2007, Kouseisha Kouseikaku.

page top

### 3-8. Marquardt method

The Newton method is a basic method for the optimization of the multivariable model

where Y is the objective function and θi is a parameter. These are non-linear equations. The Newton method to solve these equations becomes an iteration method by using the following linear equations:

where H is the n-th square matrix called a Hessian, Δθ is the modified vector of parameters, and −f is the steepest descent vector. Although the convergence area of this is narrow, the Marquardt (1963) method was invented to enlarge the convergence area as follows:

where I is the unit matrix, and λ is an adjusting factor. When λ is large, Eq. (3.66) approaches to

This is the steepest descent method and the convergence area is not large. Therefore, we use the standard scaling of parameters (Marquardt 1963; Akamine 1988b) to enlarge the convergence area and obtain

where

Thus, only the diagonal elements of H need to be multiplied by (1 + λ) for the standard scaling of parameters. When λ is large, Eq. (3.68) approaches

This is the Jacobi method (Akamine 1987b). Then the Marquardt method is a mixture of the Newton and the Jacobi methods, not a mixture of the Newton and the steepest descent methods. In the first half of the iteration, λ must be large enough to approach the Jacobi method. In the second half, λ must be small enough to approach the Newton method.

The procedure of the Marquardt method is as follows:

1. Let the initial value be λ(0) = 1, and solve Eq. (3.68).
2. When Y is decreased, let λ(t+1) = λ(t)/2, and repeat the procedure (1)–(2).
3. When Y is not decreased, let λ(t) ← 2λ(t) and solve Eq. (3.68) again. If Y is not decreased after 10 times repeating the same procedure, either the solution has been obtained or the initial values of the parameters are not adequate.

Figure 7 shows the image of the convergence (Akamine 1987a). The disadvantage point of this method is the calculation amount of the Hessian is large.

Fig. 7. Image of the convergence of the Marquardt method. M is the starting point, G is the goal, and S is the starting point of the steepest descent method which does not converge to the goal. Redrawn with permission after Akamine T. Suisan Shigen Kaiseki no Kiso, p. 11, Fig. 1.6, Kouseisha Kouseikaku, Tokyo, 2007. © 2007, Kouseisha Kouseikaku Co., Ltd.

In numerical calculation, condition (3.2) is treated by using the following methods: One way is the elimination of pn as

another way is the transformation of pi to Ki as

These ways change the conditional optimization to a normal optimization.

page top

### 4. Estimation of the population size

The population size estimation is important for fisheries stock assessment. In this section, we discuss the basic methods for estimation of the population size which are based on random sampling. These are the Petersen method, the quadrat method, and the DeLury removal method. Their statistical models are based on the binomial distribution. The classical Bayesian statistics and sum formulae of the binomial and hyper-geometric distributions are also explained. The graphic function of spread-sheet software helps to understand these methods.

page top

### 4-1. Petersen method

This is a basic method to estimate the number of fish by using the mark–recapture experiments. Let N be the total number of individuals to estimate, M be the number of marked in N, n be the total number of captured, and r be the number of mark–recaptures in n. Kuno (1986) showed the point estimation

its variance

and its 95% confidence interval

However, the precision of this interval is not high because it uses only propagation of errors.

The strict statistical model for the Petersen method is the hyper-geometric distribution

where

When N and M is larger than n and r greatly, this approximates to the binomial distribution as follows:

where p = M/N is the mark ratio.

If p in the binomial distribution model is estimated, then we can obtain the total number of individuals as N = M/p. In practice, the approximation to the normal distribution easily gives the interval estimation. According to the normal distribution, the following equation gives the confidence interval of p:

where z is the standardized variable. This is rewritten as

This is a quadratic equation of p. Akamine (2002) showed the roots of this equation with the continuity correction are

The reason that the sign in front of z is inversed is explained below. The roots without the continuity correction are

This equation is shown in Yamada and Kitada (1997). When n is large, these are approximated to

This equation is shown in many statistics textbooks.

(Example 4) Estimation of the 95% interval of N when M = 60, n = 141, and r = 11.

When z = 1.96, we can obtain the 95% confidence interval. The solutions resulting from Eq. (4.11) are N = 490.67 and 1778.01, thus the 95% confidence interval of N is [491, 1778]. The solutions resulting from Eq. (4.10) are N = 446.78 and 1360.00, thus that of N is [447, 1360]. The solutions resulting from Eq. (4.9) are N = 461.63 and 1284.13, thus that of N is [462, 1284]. The higher the precision becomes, the shorter the interval becomes. The true value of this experiment is N = 1200 and p = 0.05 (Yamada and Kitada 1997). Let r = 11.5 with the continuity correction, then we obtain z = 1.72 < 1.96. Thus, this experiment is successful and the confidence interval must include the true value. The solution from Eq. (4.9) is the best because it includes the true value and it is the shortest. When the approximation error of the binomial distribution to the normal distribution is not small, the solution from Eq. (4.10) may be safer than Eq. (4.9). On the other hand, the solutions from Eq. (4.3) are N = 374.70 and 1163.48, thus the 95% confidence interval of N is [375, 1163]. This interval does not include the true value.

The traditional interpretation of the confidence interval is as follows: The 95% confidence interval includes the true value about 95 times in 100 trials. Figure 8 shows the 95% confidence interval of the Petersen method. This is called the 95% confidence belt. Let the true value p0 be the dotted line. The data r0 is included in the 95% probability interval (the dotted line in the 95% confidence belt), and its confidence interval of p (the line segment at r = r0) includes the true value p0. On the other hand, the data r1 is not included in the 95% probability interval, and its confidence interval (the line segment at r = r1) does not include the true value p0. For interval estimation, the probability for horizontal (r) is only used, and the probability for vertical (p) is not used in traditional statistics.

The reason that the sign in front of z in Eq. (4.9) is inversed is as follows: Figure 8 shows the two curves of Eq. (4.9). The right curve gives the upper limit of r and the left one gives the lower limit. For the continuity correction, we plus 0.5 to r in the upper limit, and minus 0.5 in the lower limit. Figure 8 shows that the lower limit of p uses the upper limit of r, and the upper limit of p uses the lower limit of r. Thus, the sign in front of z in Eq. (4.9) is inversed. On the other hand, Takeuchi and Fujino (1981) used the equation in which the sign in front of z in Eq. (4.9) is not inversed, because they gave rejection regions, not a confidence interval.

Fig. 8. The 95% confidence belt of the Petersen method. Redrawn with permission after Akamine T. Suisan Shigen Kaiseki no Kiso, p. 54, Fig. 3.5, Kouseisha Kouseikaku, Tokyo, 2007. © 2007, Kouseisha Kouseikaku Co., Ltd.

page top

### 4-2. Bayesian statistical method for the Petersen method

The classical Bayesian statistics method is easily applied for this model. Bayesian statistics is based on the following theorem: The posterior distribution of θ is

where X is data, θ is parameter, L(θ|X) is likelihood, and π(θ) is the prior distribution of θ. This equation means the following proportional expression:

All methods which treat parameter probability are involved in Bayesian statistics. Equation (4.12) is also regarded as the weighted likelihood.

HDR (Highest Density Region) can be defined as follows: Let

where P is probability, x is the inner point of HDR, and y is the outside point of HDR. Thus, the area of HDR is the minimum. It is not always a continuous area for the mixed distribution. The following function can be defined as

This function gives the (1 − α) probability interval which is determined by h as

The favorable prior distribution of p is the uniform distribution from p = 0 to 1 for the Petersen method. This is the first example of Bayesian statistics in history. Let π(p) = 1. Thus, we obtain the posterior distribution

because

The following formula means that the upper probability of r and the lower probability of p are almost equal

When the binomial distribution approaches to the Poisson distribution not the normal distribution, the Bayesian result is not equal to the traditional one, because the Poisson distribution is not symmetric. When the rejection region of the normal distribution is 2.5% regions of both sides, and that of the Poisson distribution is an almost 5% region of the right side, so the results of the Bayesian method is different from the traditional one.

Solving Example 4 by using the Bayesian statistics, we obtain Fig. 9 as the posterior distribution of p with the prior distribution π(p) = 1 and the step-wise 0.005. If h = 0.025, we obtain the probability interval of p as [0.045, 0.130], and the total probability 95.32% > 95%. Removing p = 0.130 which has the minimum probability, we have the interval [0.045, 0.125] with the total probability 94.04% < 95%. In this case, the 95% probability interval of N is [480, 1333]. If p = 0.130 is included, the lower limit of N is 462, and becomes similar to the confidence interval of Eq. (4.9). When we use smaller steps, we obtain a higher precision interval.

Fig. 9. The posterior distribution of the Petersen method. Redrawn with permission after Akamine T. Suisan Shigen Kaiseki no Kiso, p. 56, Fig. 3.6, Kouseisha Kouseikaku, Tokyo, 2007. © 2007, Kouseisha Kouseikaku Co., Ltd.

page top

### 4-3. Bayesian statistical method by using the hyper-geometric distribution

The strict statistical model for the Petersen method is the hyper-geometric distribution explained above. In the binomial distribution model, we use the uniform distribution for the prior distribution of p. Therefore, the prior distribution of N is derived as follows: Let p = f(N). Thus,

Then the favorable prior distribution of N is

because the following equation holds:

Therefore, the posterior distribution of N is

The following equation means that the upper probability of r and the upper probability of N are almost equal

This model is just a simple model in which the binomial distribution is transformed to the hyper-geometric distribution. Thus, (n + 1)HG(r, n, M, N) must be the probability density in the posterior distribution. When the prior distribution is not a uniform distribution, the posterior probability is not equal to the likelihood. We had better make the probability interval of the posterior distribution be almost equal to the confidence interval of the traditional statistics (Akamine 1989b).

page top

### 4-4. Quadrat method

This is a basic method to estimate the total number of fish, plankton, or benthos by counting the number of individuals in a part (Akamine 1981, 2002). When fish distribute at random, the number of individuals caught by a part is modeled as the binomial distribution:

where n is the total number, p is the fishing ratio, and r is the catch. This is the statistical model to estimate n by using p and r. Although the point estimation is easy as n = r/p, the interval estimation is not so easy. We obtained the following estimation:

This is similar to Eq. (4.3), so the precision is not high.

In practice, the approximation to the normal distribution gives the interval estimation. Equation (4.8) is also a quadratic equation of n. The roots of Eq. (4.8) with the continuity correction are

The reason that the sign in front of z is inversed is as same as for the Petersen method. Figure 10 shows two curves of these, they are called the 95% confidence belt, when z = 1.96 and p = 0.2.

Fig. 10. The 95% confidence belt of the quadrat method. Redrawn with permission after Akamine T. Suisan Shigen Kaiseki no Kiso, p. 51, Fig. 3.3, Kouseisha Kouseikaku, Tokyo, 2007. © 2007, Kouseisha Kouseikaku Co., Ltd.

page top

### 4-5. Bayesian statistical method for the quadrat method

The favorable prior distribution of n is the uniform distribution from n = r to ∞ for the quadrat method. Let the prior distribution be π(n) = ε >0 (Mangel and Beder 1985; Akamine 1989a). Thus, we obtain the posterior distribution

because

(Example 5) Estimation of the 95% interval of n when r = 5 and p = 0.1.

The solutions of Eq. (4.27) are n = 25.38, 105.35. Thus, the 95% confidence interval of n is [26, 105].

On the other hand, the posterior distribution of n is shown in Fig. 11 by using the uniform distribution for the prior distribution in Bayesian statistics. This is the vertical probability of Fig. 10. When h = 0.0025, the interval of n is [19, 105]. The total probability of this case is 95.16%, which is larger than 95%. The total probability without n = 105 (the probability of which is minimum) is 94.91%, thus the 95% probability interval of n is [19, 104]. This interval is similar to the traditional confidence interval. The following formula supports this phenomenon.

This formula suggests that the upper probability of r and the lower probability of n are almost equal. When the binomial distribution approaches to the Poisson distribution, the Bayesian estimate and the traditional one are not equal.

Fig. 11. The posterior distribution of the quadrat method. Redrawn with permission after Akamine T. Suisan Shigen Kaiseki no Kiso, p. 52, Fig. 3.4, Kouseisha Kouseikaku, Tokyo, 2007. © 2007, Kouseisha Kouseikaku Co., Ltd.

page top

### 4-6. DeLury removal method

This method gives the estimation of the initial population size N0 and the catchability coefficient q at the same time. This method is modeled in two ways as a linear regression and a product of binomial distributions. Although Kuno (1986) recommended the former in practice, Akamine et al. (1992) recommended the expanded least-squares method based on the latter. Their objective function is given by using the approximation to the normal distribution as follows:

where Ct is the number of fish caught, Nt is the population size, pt is the removal ratio, and Xt is the fishing effort at t-th time. This is the expanded least-squares method. Therefore, the distribution of the minimum value of Y with respect to N0 and q is χ2(n − 1). Although this equation was deduced by Zippin (1963) for the minimum χ2 method, Akamine et al. (1992) obtained it from the likelihood ratio test. However, the confidence region of this model is very wide.

On the other hand, the product of binomial distributions is modified to the multinomial distribution. Thus, the other objective function is given as

The distribution of the minimum value of this with respect to N0 and q is also χ2(n− 1). Nowadays, the confidence region of N0 and q with χ2(2) is given by computer software.

(Example 6) Analysis of the data in Table 5.

The optimization method of spread-sheet software gives the minimum value of Eq. (4.31) as N0 = 1371.4, q = 0.0010651, and Ymin = 10.428. This minimum value of Y is smaller than χ2(13) = 22.36 (5%). Thus, this data passes the test of fitness. We can obtain the 95% confidence region by using the optimization method with the condition Y0 = Ymin + χ2(2) = 16.419, where χ2(2) = 5.991 (5%). This contour of Y0 is shown in Fig. 12. The maximum N0 is 2171.9 (q = 0.00060), and minimum is 1074.9 (q = 0.00150). Therefore, the confidence interval of N0 is [1075, 2171]. It is important that the upper range of the confidence interval is larger than the lower one.

Table 5. Data for the DeLury removal method. Data were excerpted with permission from Nippon Suisan Gakkaishi, 24, Nose Y. On the confidence limits corresponding to the estimate obtained by DeLury's logarithmic catch–effort method. 953–956, 1959. © 1959, the Japanese Society of Fisheries Science.

Fig. 12. The 95% confidence area of the DeLury removal method. Redrawn with permission after Akamine T. Suisan Shigen Kaiseki no Kiso, p. 61, Fig. 3.8, Kouseisha Kouseikaku, Tokyo, 2007. © 2007, Kouseisha Kouseikaku Co., Ltd.

page top

### 4-7. Proof of the sum formulae of the binomial distribution and the hyper-geometric distribution

Akamine (2002) proved some formulae as follows: The basic formula is

The sum of both sides is

The recursion of this gives Eq. (4.30). The limit is Eq. (4.29). On the other hand, integral by parts gives

The limit is Eq. (4.18), and Eqs. (4.34) and (4.35) give Eq. (4.19).

For the hyper-geometric distribution, as the same way as the binomial distribution, we obtain the following formulae:

and

Although many formulae of the hyper-geometric distribution can be derived by the use of the recurrence formula, Eqs. (4.22) and (4.24) are derived by the use of the summation by parts. The difference is defined as follows:

Hence, we obtain

Let Δf(x) = g(x), this gives the indefinite summation

where Δ−1 is the summation operator, c is the summation constant. Therefore, the definite summation is denoted as

The difference of the product is

Thus, we obtain the following two formulae of the summation by parts:

The differences of the combination with respect to x are

and

The inverse operation of this is

Formula (4.22) can be modified as

We use formula (4.45) for S.

Thus, formula (4.22) is proved. Next, the left side of formula (4.24) can be modified as

We use formula (4.46) for S*.

Therefore, we obtain the following equation:

Finally, Eqs. (4.54), (4.37) and (4.38) give Eq. (4.24).

page top

### 5. Survival models

The virtual population analysis (VPA) is a typical survival model in fish population dynamics (Megrey 1989). Most of them use mortality coefficients. On the other hand, the Leslie matrix model is a famous management model in agriculture. The VPA model using mortality rates is useful because it is related to the Leslie matrix model, and linear programming. In this section, these models are presented as simple styles, because simple ones are most useful for many applications.

page top

### 5-1. VPA

The standard fishing models are

where N(t) is the number of fish, C(t) is the number of fish caught, Z is the total mortality coefficient, F is the fishing mortality coefficient, and M is the natural mortality coefficient. The resulting equation from these two equations is

Aksland (1994) and Hiramatsu (1995) obtained the analytical solution of this as

On the other hand, let U(t) be the number of natural mortality and we have

Hence, Eqs. (5.1), (5.2) and (5.5) give a simple relationship

Equations (5.1) and (5.2) also give the relationship

When F and M are constants, the integrals of Eqs. (5.1) and (5.2) are

where N0 = N(0). VPA is a simple estimating method of N from C by using these fishing equations.

The definite calculation of VPA is as follows: The constant value of M is known. The basic equations are

where Ni = N(i), and

Substituting Eq. (5.10) into Eq. (5.11) gives

At first, we solve this non-linear equation for Fi, and then we obtain Ni of Eq. (5.10). The recurrence of these procedures gives all Ni. This method needs the number of fish in the oldest class NT which is called the terminal N. In practice, we use the terminal F instead of N. The following Ishioka and Kishida (1985) iteration method for Eq. (5.13) has a high precision and converges at about three times recurrence. Equation (5.13) can be modified as follows:

This is a special case of Eq. (5.4). The Ishioka and Kishida iteration method is

On the other hand, a few approximations are used for this non-linear equation. The famous Pope method is

Thus, Eq. (5.14) is rewritten as

This is an impulse fishing model at the middle point in the fishing period. The fishing equation of this model is

The following Nagai (1980) method has a higher precision than the Pope method. The Nagai method is

This is also an impulse fishing model in which the fishing point is a little before that of the Pope method. This approximation is also regarded as

where

The last function is the generating function of Bernoulli numbers. The complete solution of Eq. (5.20) is only the exponential function.

Akamine (2001) defined the impulse fishing as follows: Let

Hence, we have

When ba, we obtain

where δ is the Dirac delta function, and

Equation (5.18) is a special case of this.

page top

### 5-2. VPA using mortality rates

The following simple model is used for whales:

where E is the fishing mortality rate. This is an impulse fishing model at the first point in the fishing period. The model becomes the same as the Pope method by sliding a fishing time for a half period.

The discrete fishing models are

where D is the natural mortality rate. Thus the fishing models using mortality rates instead of mortality coefficients are as follows (Akamine 1995a, 1995d):

where s is the survival rate. Hence, forward-calculation is

and back-calculation is

This is rewritten as

where vi = 1/(1 − Di). This is equivalent to Eq. (5.4) or Eq. (5.14) essentially. When D is constant, we have

The relationships of mortality rates and mortality coefficients are

The inverse transformations are

In these transformations, the following relationships hold:

These are shown in Fig. 13, in which solid and broken lines have a one-to-one correspondence.

Fig. 13. The transformation of the coefficients and rates. Redrawn after Bull. Jpn. Soc. Fish. Oceanogr., 59, Fig. 1, Akamine T. Introduction to cohort analysis (VPA), 424–437, 1995, with permission from the Japanese Society of Fisheries Oceanography.

page top

### 5-3. Leslie matrix model

This is a famous discrete survival model represented as

The representation by the components of this is

where Ni(j) is the number of i-th year old fish in j year, Ri is the reproduction rate and si is the survival rate of i-th year old fish. The resolution of this Leslie matrix is

where ri is the number of eggs reproduced by one individual and Ri = risi. Thus, the forward-calculation of VPA is represented by matrix as follows (Akamine 1995b):

page top

### 5-4. Linear programming for fishing equations

Akamine (1996, 1997) showed the optimization by linear programming for fishing. He rewrote Eq. (5.32) as

where mi is constant. Resulting from Eq. (5.30) is

Definition of the reproduction is

and that of the yield is

Let the fishing mortality rates be control variables as

Then we obtain the inequalities

The objective function (5.49) is rewritten as

Therefore, we can use the linear programming for the objective function (5.52) under conditions (5.48) and (5.51).

page top

### 6. Summary

In this monograph, the author has presented four subjects. The first is the new standard growth formula and its statistical methods for estimation and test. The second subject is the body-size composition analysis in which we use the maximum likelihood method for parameter estimation on a mixture of normal distributions. The third subject is the estimation of population size in which the Petersen method, the quadrat method, and the DeLury removal method are discussed. The last subject is the survival models in which VPA, the Leslie matrix model, and the linear programming are explained. We can use all the methods in this monograph by using spread-sheet software.

In Section 2, a new standard growth formula based on the Richards growth formula and a periodic function

is presented. Parameter estimation and test for this formula is based on the weighted least-squares method. The Ohnishi and Akamine growth formula, which is an implicit function, and the Awaya method to estimate parameters of implicit models are introduced. The generalized reproduction model which is similar to the Richards growth formula is also explained.

In Section 3, the Hasselblad method for estimation on a mixture of normal distributions is discussed. This algorithm is obtained by using the EM algorithm naturally, and the convergence of its objective function is proved by using the K–L information quantity. The sufficient condition for convergence of parameters is as follows:

where L is the Lipschitz constant. On the other hand, the modification of the Jacobi method is

This shows the convergence speed of the Hasselblad method, which is the numerator of the right side, is similar to the Jacobi method.

In Section 4, the traditional methods in which the binomial distribution is approximated to the normal distribution for the Petersen method and the quadrat method are explained. The Bayesian statistical methods for these are also presented. The hyper-geometric distribution is used for the Petersen method in Bayesian, and its sum formulae are proved by using the formulae of summation by parts. The following formulae show the relationship between the traditional method and the Bayesian method:

and

The DeLury removal method whose objective function is defined as the expanded weighted least-squares

is also explained.

In Section 5, the Pope, Nagai methods, and Ishioka and Kishida iteration method for VPA are explained. The simple VPA method using mortality rates instead of mortality coefficients is presented as follows:

The transformation between mortality rates and mortality coefficients are

The Leslie matrix model related to VPA, and the linear programming using mortality rates are also presented.

page top

### Acknowledgements

The author thanks the late Mr. Fumihiko Kato and the late Mr. Kiyohide Ishioka for their encouragement. The author also thanks the two anonymous referees for their useful comments on earlier drafts of the manuscript.

page top

### References

Aizawa Y, Takiguchi N. Consideration of the methods for estimating the age-composition from length frequency data with MS-Excel. Bull. Jpn. Soc. Fish. Oceanogr. 1999; 63: 205–214 (in Japanese).

Akamine T. A handy method for selecting and counting bivalve larvae at planktonic stage. Bull. Japan Sea Natl. Fish. Res. Inst. 1981; 32: 77–81 (in Japanese).

Akamine T. A BASIC program to analyse the polymodal frequency distribution into normal distributions. Bull. Japan Sea Natl. Fish. Res. Inst. 1982; 33: 163–166 (in Japanese).

Akamine T. The BASIC program to analyse the polymodal frequency distribution into normal distributions with Marquardt's method. Bull. Japan Sea Natl. Fish. Res. Inst. 1984; 34: 53–60 (in Japanese).

Akamine T. Consideration of the BASIC programs to analyse the polymodal frequency distribution into normal distributions. Bull. Japan Sea Natl. Fish. Res. Inst. 1985; 35: 129–159 (in Japanese).

Akamine T. Expansion of growth curves using a periodic function and BASIC programs by Marquardt's method. Bull. Japan Sea Natl. Fish. Res. Inst. 1986; 36: 77–107.

Akamine T. A solution of the multi-cohort model by Marquardt's method. Bull. Japan Sea Natl. Fish. Res. Inst. 1987a; 37: 225–257.

Akamine T. Comparison of algorithms of several methods for estimating parameters of a mixture of normal distributions. Bull. Japan Sea Natl. Fish. Res. Inst. 1987b; 37: 259–277.

Akamine T. Evaluation of error caused by histogram on estimation of parameters for a mixture of normal distributions. Bull. Japan Sea Natl. Fish. Res. Inst. 1988a; 38: 171–185.

Akamine T. Estimation of parameters for Richards model. Bull. Japan Sea Natl. Fish. Res. Inst. 1988b; 38: 187–200.

page top

Akamine T. An interval estimation for extraction using Bayesian statistics. Bull. Japan Sea Natl. Fish. Res. Inst. 1989a; 39: 9–17.

Akamine T. An interval estimation for the Petersen method using Bayesian statistics. Bull. Japan Sea Natl. Fish. Res. Inst. 1989b; 39: 19–35.

Akamine T. A new standard formula for seasonal growth of fish in population dynamics. Nippon Suisan Gakkaishi 1993; 59: 1857–1863.

Akamine T. Definition of the exponential function and its application for population dynamics. Bull. Jpn. Soc. Fish. Oceanogr. 1994; 58: 317–320 (in Japanese).

Akamine T. Mathematical approach to virtual population analysis. Fish. Sci. 1995a; 61: 167.

Akamine T. Relationship between Leslie matrix and cohort analysis. Fish. Sci. 1995b; 61: 722.

Akamine T. A mathematical study of fish growth formula in population dynamics. Bull. Natl. Res. Inst. Fish. Sci. 1995c; 7: 189–263 (in Japanese).

Akamine T. Introduction to cohort analysis (VPA). Bull. Jpn. Soc. Fish. Oceanogr. 1995d; 59: 424–437 (in Japanese).

Akamine T. Discrete fishing equations and linear programming. Bull. Jpn. Soc. Fish. Oceanogr. 1996; 60: 252–258 (in Japanese).

Akamine T. Optimum control theory for the dynamic pool model. Bull. Natl. Res. Inst. Fish. Sci. 1997; 10: 135–167 (in Japanese).

page top

Akamine T. Convergence of Hasselblad method for estimating proportions of a mixture of probability distributions. Fish. Sci. 1998; 64: 997–998.

Akamine T. Convergence of Hasselblad method for estimating parameters of a mixture of normal distributions. Bull. Natl. Res. Inst. Fish. Sci. 1999; 14: 49–58 (in Japanese).

Akamine T. Mathematical study for approximations and iterations of the virtual population analysis. Bull. Natl. Res. Inst. Fish. Sci. 2001; 16: 1–16 (in Japanese).

Akamine T. Comparison between conventional and Bayesian statistics in interval estimation for quadrat method and Petersen method. Bull. Fish. Res. Agen. 2002; 2: 25–34 (in Japanese).

Akamine T. Statistical test and model selection of fish growth formula. Bull. Jpn. Soc. Fish. Oceanogr. 2004; 68: 44–51 (in Japanese).

Akamine T. Mixture of normal distributions and EM algorithm. Bull. Jpn. Soc. Fish. Oceanogr. 2005; 69: 174–183 (in Japanese).

Akamine T. Suisan Sigen Kaiseki no Kiso. Kouseisha Kouseikaku, Tokyo. 2007; 115 pp.

Akamine T, Matsumiya Y. Mathematical analysis of age–length key method for estimating age composition from length composition. Bull. Japan Sea Natl. Fish. Res. Inst. 1992; 42: 17–24.

Akamine T, Kishino H, Hiramatsu K. Non-biased interval estimation of Leslie's removal method. Bull. Japan Sea Natl. Fish. Res. Inst. 1992; 42: 25–39.

Aksland M. A general cohort analysis method. Biometrics 1994; 50: 917–932.

page top

Appeldoorn RS. Modification of a seasonally oscillating growth function for use with mark–recapture data. J. Cos. Int. Explor. Mer 1987; 43: 194–198.

Awaya T. Two-dimensional curve fitting in counting experiments. Nucl. Instr. Meth. 1983; 212: 311–317.

Awaya T. Data Kaiseki. Gakkai Shuppan Center, Tokyo. 1991; 270 pp.

Cloern JE, Nichols FH. A von Bertalanffy growth model with a seasonally varying coefficient. J. Fish. Res. Board Can. 1978; 35: 1479–1482.

Dark TA. Age and growth of Pacific hake, Merluccius productus. Fish. Bull. 1975; 73: 336–355.

France J, Thornley JHM. Mathematical Models in Agriculture. Butterworth, London. 1984; pp. 75–94.

Gorie S. Estimation of the parameters in von Bertalanffy's growth formula by MS-Excel. Suisanzoshoku 2001; 49: 519–527 (in Japanese).

Gorie S. Estimation of parameters in a mixture of normal distributions from length frequency composition and growth formula by MS-Excel. Suisanzoshoku 2002; 50: 243–249 (in Japanese).

Group of Statistics, Pelagic Resource Section, SEIKAI Regional Fisheries Research Laboratory. Data of the investigation of September to November in 1950. Report of the Investigations of Trawl Fishery in the China Sea, No. 2. SEIKAI Regional Fisheries Research Laboratory, Nagasaki, Japan. 1951; pp. 1–28.

Haddon M. Modelling and Quantitative Methods in Fisheries. Chapman and Hall, New York. 2001; 406 pp.

page top

Hasselblad V. Estimation of parameters for a mixture of normal distributions. Technometrics 1966; 8: 431–444.

Hiramatsu K. A theoretical study of equations used in virtual population analysis. Fish. Sci. 1995; 61: 752–754.

Iri M. Suuchi Keisan. Asakura Publishing, Tokyo. 1981; 173 pp.

Ishioka K, Kishida T. Studies on the algorithm for solving the catch equation and its accuracy. Bull. Nansei Reg. Fish. Res. Lab. 1985; 19: 111–120 (in Japanese).

Kamijoh Y, Yokomatsu Y, Andoh K. Hamaguri no sigen baiyou gijutu kaihatu kennkyuu houkokusho. Oita Prefectural Agriculture, Forestry and Fisheries Research Center. 1985; 54 pp.

Kimura DK. Likelihood methods for the von Bertalanffy growth curve. Fish. Bull. 1980; 77: 765–776.

Kiso K, Akamine T, Ohnishi S, Matsumiya Y. Mathematical examinations of the growth of sea-run and fluviatile forms of female Masu salmon Oncorhynchus masou in rivers of the southern Sanriku district, Honshu, Japan. Nippon Suisan Gakkaishi 1992; 58: 1779–1784.

Kuno E. Doubutsu no Kotaigun Doutai Kenkyuu Hou 1, Kotaisuu Suitei Hou. Kyoritsu Shuppan, Tokyo. 1986; 114 pp.

Mangel M, Beder JH. Search and stock depletion: theory and applications. Can. J. Fish. Aquat. Sci. 1985; 42: 150–163.

Marquardt DW. An algorithm for least-squares estimation of nonlinear parameters. J. Soc. Indust. Appl. Math. 1963; 11: 431–441.

page top

Megrey BA. Review and comparison of age-structured stock assessment models from theoretical and applied points of view. American Fisheries Society Symposium 1989; 6: 8–48.

Murata N, Ikeda S. Neural network to EM. In: Watanabe M, Yamaguchi K (eds). EM Algorithm to Fukanzen Data no Shomondai. Taga Shuppan, Tokyo. 2000; pp. 155–188.

Nagai T. Cohort analysis nyuumon. Enyou Suiken News 1980; 38: 10–12 (in Japanese).

Nose Y. On the confidence limits corresponding to the estimate obtained by DeLury's logarithmic catch–effort method. Nippon Suisan Gakkaishi 1959; 24: 953–956 (in Japanese).

Ohnishi S, Akamine T. Extension of von Bertalanffy growth model incorporating growth patterns of soft and hard tissues in bivalve molluscs. Fish. Sci. 2006; 72: 787–795.

Pauly D, David N. ELEFAN 1, a BASIC program for the objective extraction of growth parameters from length-frequency data. Meeresforsch. 1981; 28: 205–211.

Pauly D, Gaschutz G. A simple method for fitting oscillating length growth data with a program for pocket calculators. I.C.E.S. CM 1979/G/24, Demersal Fish Committee. 1979; 26 pp.

Pauly D, Soriano-Bartz M, Moreau J, Jarre-Teichmann A. A new model accounting for seasonal cessation of growth in fishes. Aust. J. Mar. Freshwater Res. 1992; 43: 1151–1156.

Pitcher TJ, MacDonald PDM. Two models for seasonal growth in fishes. J. Appl. Ecol. 1973; 10: 599–606.

Richards FJ. A flexible growth function for empirical use. J. Exp. Bot. 1959; 10: 290–300.

page top

Sakamoto Y, Ishiguro M, Kitagawa G. Akaike Information Criterion Statistics. D.Reidel Publishing, Tokyo. 1986; xix+290 pp.

Schnute J. A versatile growth model with statistically stable parameters. Can. J. Fish. Aquat. Sci. 1981; 38: 1128–1140.

Schnute J. A general theory for analysis of catch and effort data. Can. J. Fish. Aquat. Sci. 1985; 42: 414–429.

Somers IF. On a seasonally oscillating growth function. Fishbite 1988; 6: 8–11.

Takeuchi K, Fujino K. 2-Kou Bunpu to Poisson Bunpu. University of Tokyo Press, Tokyo. 1981; 262 pp.

Tanaka S. A method of analyzing the polymodal frequency distribution and its application to the length distribution of porgy. Bull. Tokai Reg. Fish. Res. Lab. 1956; 14: 1–13 (in Japanese).

Taylor CC. Growth equations with metabolic parameters. J. Cons. Int. Explor. Mer 1962; 27: 270–286.

Yamada S, Kitada S. Seibutsu Sigen Toukeigaku. Seizando, Tokyo. 1997; 263 pp.

Zippin C. An evaluation of the removal method of estimating animal populations. Biometrics 1963; 12: 163–169.

page top

### List of Figures

Table 1. The growth data of the clam (Meretrix lusoria). SL is the Shell length (mm). Data were excerpted from Kamijoh et al. Hamaguri no shigen baiyou gijutu kaihatu kenkyu houkokusho, 1985, with permission from Fisheries Research Institute, Oita Prefectural Agriculture, Forestry and Fisheries Research Center.

Table 2. The body-length data (cm) of the Pacific hake (Merluccius productus). Data were excerpted with permission from Fishery Bulletin, 73, Dark TA. Age and growth of Pacific hake, Merluccius productus. 336–355, 1975.

Table 3. The body-length composition data of the porgy (Dentex tumifrons). BL is the body length (cm). Data were excerpted with kind permission of Dr. Shoichi Tanaka from Table 2 in Tanaka (1956) with its original data source of Group of Statistics, Pelagic Resource Section, SEIKAI Regional Fisheries Research Laboratory, Data of the investigation of September to November in 1950, Report of the Investigations of Trawl Fishery in the China Sea, No. 2, pp. 1–28, 1951. SEIKAI Regional Fisheries Research Laboratory, Nagasaki, Japan.

Table 4. The result of the estimation for the mixture of normal distributions. The range of the body length is 7–32 cm. Data were excerpted with permission from Akamine T. Suisan Shigen Kaiseki no Kiso, Table 1.2, Kouseisha Kouseikaku, Tokyo, 2007. © 2007, Kouseisha Kouseikaku.

Table 5. Data for the DeLury removal method. Data were excerpted with permission from Nippon Suisan Gakkaishi, 24, Nose Y. On the confidence limits corresponding to the estimate obtained by DeLury's logarithmic catch–effort method. 953–956, 1959. © 1959, the Japanese Society of Fisheries Science.

page top

Fig. 1. Graphs of the Richards growth formulae. Figures are values of r. Revised from Bull. Japan Sea Natl. Fish. Res. Inst., 38, Fig. 2, Akamine T. Estimation of parameters for Richards model. 187–200, 1988, with permission from Japan Sea National Fisheries Research Institute, Fisheries Research Agency.

Fig. 2. The growth curve of the clam. Black circles are SL, and squares are S.D. Revised with permission from Akamine T. Taicho-sosei data karano nenrei to seicho no suitei. In: Shigen Hyoka notameno Suuchi Kaiseki, Suisangaku Series 66, p. 127, Fig. 8.6, Kouseisha Kouseikaku, Tokyo, 1987. © 1987, the Japanese Society of Fisheries Science.

Fig. 3. Black circles and the thick line are female data and white circles and the thin line are male. Redrawn with permission after Fishery Bulletin, 77, Fig. 1, Kimura DK. Likelihood methods for the von Bertalanffy growth curve. 765–775, 1980.

Fig. 4. The 95% confidence interval of the growth curve for female Pacific hake. Two thick lines show growth curves for the minimum and maximum k. Redrawn with permission after Akamine T, Suisan Shigen Kaiseki no Kiso, p. 32, Fig. 2.6. Kouseisha Kouseikaku, Tokyo, 2007. © 2007, Kouseisha Kouseikaku Co., Ltd.

Fig. 5. The graphs of the generalized reproduction model. Figures are values of r. Redrawn with permission after Schnute J. A general theory for analysis of catch and effort data. Can. J. Fish. Aquat. Sci., 42(3), 419–429, Fig. 2, 1985. © 2008, NRC Canada.

page top

Fig. 6. Result of the estimation (\optS{Tabe 4}). Black circles are data, thick line is a mixture, and thin lines are normal distributions. Redrawn with permission after Akamine T. Suisan Shigen Kaiseki no Kiso, p. 6, Fig. 1.3, Kouseisha Kouseikaku, Tokyo, 2007. © 2007, Kouseisha Kouseikaku Co., Ltd.

Fig. 7. Image of the convergence of the Marquardt method. M is the starting point, G is the goal, and S is the starting point of the steepest descent method which does not converge to the goal. Redrawn with permission after Akamine T. Suisan Shigen Kaiseki no Kiso, p. 11, Fig. 1.6, Kouseisha Kouseikaku, Tokyo, 2007. © 2007, Kouseisha Kouseikaku Co., Ltd.

Fig. 8. The 95% confidence belt of the Petersen method. Redrawn with permission after Akamine T. Suisan Shigen Kaiseki no Kiso, p. 54, Fig. 3.5, Kouseisha Kouseikaku, Tokyo, 2007. © 2007, Kouseisha Kouseikaku Co., Ltd.

Fig. 9. The posterior distribution of the Petersen method. Redrawn with permission after Akamine T. Suisan Shigen Kaiseki no Kiso, p. 56, Fig. 3.6, Kouseisha Kouseikaku, Tokyo, 2007. © 2007, Kouseisha Kouseikaku Co., Ltd.

Fig. 10. The 95% confidence belt of the quadrat method. Redrawn with permission after Akamine T. Suisan Shigen Kaiseki no Kiso, p. 51, Fig. 3.3, Kouseisha Kouseikaku, Tokyo, 2007. © 2007, Kouseisha Kouseikaku Co., Ltd.

page top

Fig. 11. The posterior distribution of the quadrat method. Redrawn with permission after Akamine T. Suisan Shigen Kaiseki no Kiso, p. 52, Fig. 3.4, Kouseisha Kouseikaku, Tokyo, 2007. © 2007, Kouseisha Kouseikaku Co., Ltd.

Fig. 12. The 95% confidence area of the DeLury removal method. Redrawn with permission after Akamine T. Suisan Shigen Kaiseki no Kiso, p. 61, Fig. 3.8, Kouseisha Kouseikaku, Tokyo, 2007. © 2007, Kouseisha Kouseikaku Co., Ltd.

Fig. 13. The transformation of the coefficients and rates. Redrawn after Bull. Jpn. Soc. Fish. Oceanogr., 59, Fig. 1, Akamine T. Introduction to cohort analysis (VPA), 424–437, 1995, with permission from the Japanese Society of Fisheries Oceanography.

page top