doi:10.5047/absm.2009.00203.0001
© 2009 TERRAPUB, Tokyo. All rights reserved.
www.terrapub.co.jp/onlinemonographs/absm/
National Research Institute of Fisheries Science, Fisheries Research Agency,
Fukuura 2-12-4, Yokohama 236-8648, Japan
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.
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
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.
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.
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 = ∞, t_{0} 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, w_{0}). 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 w_{f} = 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.
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.
Pitcher and MacDonald (1973) presented the growth formulae (PM1)
and (PM2)
where t_{r} is real time (weeks), t_{g} is growth time, t_{s} = t_{r} − s, s is the starting point for cosine, sw is the threshold value for cosine, C is the magnitude constant, and s_{1} 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 s_{2} is the time lag between temperature and growth switch. Using this function, they rewrote the growth formulae as follows (PM1):
and (PM2):
where sw_{1} is the threshold temperature, and s_{3} 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(t_{0}) ≠ 0.
Cloern and Nichols (1978) defined the differential equation
where L'_{max} is maximal growth rate when body size is minimal, and L_{max} 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 L_{min} 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.
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 a ≤ f(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 t → F(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.
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 t_{0} in Eq. (2.29) as follows:
Pauly and Gashutz (1979) also used the multiple regression to obtain values of C, k, t_{0} and t_{s} 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 t_{i}, let n_{i} be the number of specimens, be the mean, and σ^{2}_{i} be the variance. When t = t_{i} the residual sum of squares is
where the subscript i is omitted, w_{j} is the body-size of the j-th individual at time t_{i}, 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 Y_{0} is the objective function of the true values and Y_{min} 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 L_{max} 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 S_{0} is the residual squares of the true values and S_{min} 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 m_{M}, that of female data be m_{F}. 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 S_{0}. The null hypothesis must be
In this case p = 3, m = 13, F(3,10) = 3.708 (5%), S_{min} = S_{F} = 28.8, and we obtain S_{0} = 60.9 for the condition. By using the optimization method of spread-sheet software, we obtain the minimum k = 0.200 (l_{∞} = 66.3, t_{0} = −0.761), and the maximum k = 0.413 (l_{∞} = 58.3, t_{0} = 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 t_{0}, negative correlations between k and l_{∞}, and l_{∞} and t_{0} in VBGF1. The minimum and maximum estimates of l_{∞} and t_{0} 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.
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:
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), α = e^{−t} and β = b(1 − e^{−rt})/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).
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 (s_{i}, r_{i}, δ_{si}^{2}, δ_{ri}^{2}):
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 (t_{i}, μ_{i},δ_{i}^{2}) 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 v_{i} = 1/δ_{i}^{2}, and λ_{i} is the undetermined multiplier. Thus, the following simultaneous equations can be obtained:
where f_{i} = f(t_{i}, w_{i}, θ), 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 f_{i0} = f(t_{i}, w_{i0}, θ_{0}). Thus, the following equations can be obtained:
On the other hand, substituting Eq. (2.97) into Eq. (2.94) with f_{i} = f_{i0} gives
where d_{i0} = μ_{i} − w_{i0}. Substituting Eq. (2.100) into Eq. (2.101) gives
Substituting these into Eq. (2.95) with f_{i} = f_{i0} 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 w_{i}, 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 Δw_{i} from Eq. (2.100) and new values of f_{i} from Eq. (2.99). Therefore, solving Eq. (2.103) again and repeating these procedures gives the solution. The linear approximation errors of f_{i} 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 w_{i}.
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
Y_{2} is better than Y_{1} because Y_{2} 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 = ξ l^{2}, w = ηl^{3}. 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_{∞}, l_{0} = l(0). This is an implicit model, so the Awaya method is useful for parameter estimation.
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.
where x is the body-size, subscript i is the number of each distribution, and p_{i} > 0 is a mixture ratio which has the condition
Function g_{i}(x) ≥ 0 is a normal distribution:
The estimates of p_{i}, μ_{i} and σ_{i}^{2} 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 x_{1}, ..., x_{m}, 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.
Hasselblad (1966) defined the mixture ratio as follows:
In this case, the condition of the solution
is rewritten as
where
By multiplying p_{i} 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).
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 p_{i} to both sides of this equation, we obtain Eq. (3.16) and iteration (3.17).
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 z_{i}(x) as the relative frequency of age i and size x, and q_{i}(x) as the probability of age i in the size x as follows:
In E-step, z_{i}(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 p_{i}, 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 σ_{i}^{2}, 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 x ≤ x − 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 z_{i}(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 z_{i}(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.
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 i ≠ k 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
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 p_{i} in the Jacobi method is
This shows that the Hasselblad method is regarded as an approximation of the Jacobi method when Ψ_{i} is small.
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.
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:
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 p_{n} as
another way is the transformation of p_{i} to K_{i} as
These ways change the conditional optimization to a normal optimization.
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.
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 p_{0} be the dotted line. The data r_{0} 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 = r_{0}) includes the true value p_{0}. On the other hand, the data r_{1} is not included in the 95% probability interval, and its confidence interval (the line segment at r = r_{1}) does not include the true value p_{0}. 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.
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.
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).
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.
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.
This method gives the estimation of the initial population size N_{0} 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 C_{t} is the number of fish caught, N_{t} is the population size, p_{t} is the removal ratio, and X_{t} 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 N_{0} 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 N_{0} and q is also χ^{2}(n− 1). Nowadays, the confidence region of N_{0} 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 N_{0} = 1371.4, q = 0.0010651, and Y_{min} = 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 Y_{0} = Y_{min} + χ^{2}(2) = 16.419, where χ^{2}(2) = 5.991 (5%). This contour of Y_{0} is shown in Fig. 12. The maximum N_{0} is 2171.9 (q = 0.00060), and minimum is 1074.9 (q = 0.00150). Therefore, the confidence interval of N_{0} 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.
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).
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.
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 N_{0} = 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 N_{i} = N(i), and
Substituting Eq. (5.10) into Eq. (5.11) gives
At first, we solve this non-linear equation for F_{i}, and then we obtain N_{i} of Eq. (5.10). The recurrence of these procedures gives all N_{i}. This method needs the number of fish in the oldest class N_{T} 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 b → a, we obtain
where δ is the Dirac delta function, and
Equation (5.18) is a special case of this.
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 v_{i} = 1/(1 − D_{i}). 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.
The representation by the components of this is
where N_{i}(j) is the number of i-th year old fish in j year, R_{i} is the reproduction rate and s_{i} is the survival rate of i-th year old fish. The resolution of this Leslie matrix is
where r_{i} is the number of eggs reproduced by one individual and R_{i} = r_{i}s_{i}. Thus, the forward-calculation of VPA is represented by matrix as follows (Akamine 1995b):
where m_{i} 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).
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.
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.
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).
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.
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.
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.
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.
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.
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.
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.
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.
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.