Stand Diameter Distribution Modelling and Prediction Based on Richards Function

The objective of this study was to introduce application of the Richards equation on modelling and prediction of stand diameter distribution. The long-term repeated measurement data sets, consisted of 309 diameter frequency distributions from Chinese fir (Cunninghamia lanceolata) plantations in the southern China, were used. Also, 150 stands were used as fitting data, the other 159 stands were used for testing. Nonlinear regression method (NRM) or maximum likelihood estimates method (MLEM) were applied to estimate the parameters of models, and the parameter prediction method (PPM) and parameter recovery method (PRM) were used to predict the diameter distributions of unknown stands. Four main conclusions were obtained: (1) R distribution presented a more accurate simulation than three-parametric Weibull function; (2) the parameters p, q and r of R distribution proved to be its scale, location and shape parameters, and have a deep relationship with stand characteristics, which means the parameters of R distribution have good theoretical interpretation; (3) the ordinate of inflection point of R distribution has significant relativity with its skewness and kurtosis, and the fitted main distribution range for the cumulative diameter distribution of Chinese fir plantations was 0.4∼0.6; (4) the goodness-of-fit test showed diameter distributions of unknown stands can be well estimated by applying R distribution based on PRM or the combination of PPM and PRM under the condition that only quadratic mean DBH or plus stand age are known, and the non-rejection rates were near 80%, which are higher than the 72.33% non-rejection rate of three-parametric Weibull function based on the combination of PPM and PRM.


Introduction
Diameter distributions are well known and widely used for describing forest stand diameter structure [1,2]. Accurate quantification of tree characteristics permits study of the interaction among physical and physiological processes and growth. Quantification of diameter distributions over time allows the manager to relate the parameters of the distribution to stand age or stand density [3]. The stand volume characteristics are calculated using diameter distribution and tree height and volume models [4]. Growth and yield prediction based on the diameter distribution approach has also been widely used [5,6].
The methods are used to describe diameter distribution can be classified into parametric and nonparametric methods. The abovementioned functions and equations all belong to parametric methods. Nonparametric methods, like percentile prediction method [21,22] and k-nearest neighbor estimation method [23,24], do not rely on any predefined functional form and adapt to description of multimodal distributions. The disadvantage of nonparametric methods is that their high amount of required reference material is difficult to acquire and time-consuming [25].
Doubtlessly, in parametric methods, the Weibull is the most commonly used probability density function for fitting tree diameter distributions. Since three-parametric Weibull function had been derived by Weibull [26], due to the relative simplicity of expression formula and its flexibility in fitting a variety of shapes and degrees of skew, this function has proved to be a good distribution model. For the estimation of Weibull parameters, many different methods have been applied, such as moment method, maximum likelihood method, percentiles method and nonlinear regression method [5,11,27,28]. With the presentation of advanced analysis software, nonlinear regression method shows its superiority. The best advantage of this method is that parameters can simultaneously accurately be estimated. Additionally, we should acknowledge the facts that iterated function of the three-parametric Weibull distribution is not easy to converge while using Weibull to fit the diameter distribution data, and the correlativity between the parameters estimates and the whole stand characteristics become weak [29,30]. Thus we wish to explore a new diameter distribution model that overcomes the disadvantages of Weibull and has the advantages of Weibull such as theoretical meanings of parameters and high simulation precision.
The objectives of this study were to explore the application of Richards equation on modelling and prediction of stand diameter distribution, test the theoretical meanings of its parameters, and compare the properties of modelling and prediction for stands diameter distribution between Richards equation and threeparametric Weibull equation using the long-term repeated measurement data sets from Chinese fir (Cunninghamia lanceolata) plantations in southern China.

Data Source
Chinese fir (Cunninghamia lanceolata) stands located in Fenyi city, Jiangxi Province, China, experience a subtropical climate. The longitude is 114u339E, latitude 27u349N. Mean annual temperature, rainfall and evaporation are 16.8uC, 1656 mm, and 1503 mm, respectively. Chinese fir stands mentioned as follows in the location all are built and authorized by Research Institute of Forestry of Chinese Academy of Forestry and the data originated from our continuous survey. So no specific permits were required for the described field studies, and the field studies did not involve endangered or protected species.
The unthinned stands of Chinese fir were established in 1981. Planting density was limited within an optimum range according to managerial purposes. The series of stand planting densities was 1667, 3333, 5000, 6667, 10000 stems ha 21 . Every planting density had 3 designed replications. Each plot area was 0.06 ha and two adjacent plots were separated by buffer zone. All trees in each plot were marked for continuous measurement. Stem diameter at breast height (Dbh) was measured after tree height reached 1.3 m. All stands were measured every year before reaching 10 years old, and every two years after reaching 10 years old; all stands were measured 10 times. Self-thinning occurred in all stands during the experimental period. The basic information of 150 unthinned stands is described in Table 1. The database obtained from these unthinned stands was used to build models of diameter distribution.
Another database comprised of 159 diameter frequency distributions was used to test the models. The data came from a thinning study of Chinese fir plantations established in the same environment as the above-mentioned unthinned stands. Among the 159 diameter distributions, 63 distributions came from unthinned stands, the remaining came from thinned stands. Because of thinning only being viewed as a management measure that influences stands diameter distribution as same as density, the thinned stands were included in the test database. All thinning was from below. Each plot area was 0.05 ha. The basic information of used stands data for model evaluation is described in Table 2.

Computation of Observed Cumulative Diameter Distribution
The diameter classes applied were 2 cm wide. Diameter class, k, is defined in absolute scale (e.g., 1-2.9 cm for k = 2, 3-4.9 for k = 4, etc.), namely, diameter class k is the midpoint value of the absolute scale. The frequency of stems in diameter class k at stand i is given by: where N ki is the number of trees of diameter class k at stand i (i = 1, 2, …, 159), and N i is the total number of trees in stand i. The cumulative frequency of stems in diameter class k at stand i can be obtained by: C ki~F1i zF 2i z ::: zF (k{1)i zF ki Table 3 shows the basic mathematical characteristics of Richards and Weibull equations. In the expression of Weibull equations, a is the location parameter, b is the scale parameter, and c is the shape parameter. What, then, has caused the iterated function to not easily converge while employing SAS's nonlinear regression method [31] to estimate the parameters of Weibull and, subsequently, for the correlativity between the parameter estimates and the whole stand characteristics to be weak? When the expression formula of the Weibull equation is analyzed, it is found the equation is meaningful only when (x{a)=b §0, because there is b.0, so x §a, this means the value of parameter a should be smaller than the midpoint of minimum diameter class. If xva, Weibull equation is meaningless. Furthermore, it is difficult to give a suitable initial estimation for the location parameter a [32], and the ultimate estimation of the location parameter cannot exceed the two neighboring diameter class of any given initial estimation, which can greatly increase the times of iteration. Obviously, it is the location of shape parameter which leads to the above-mentioned shortcomings of Weibull. Richards equation have been found to be suitable to fitting the sigmoid data sets [15]. In the beginning, the curve is concave up, while in later life it becomes convex [33,34]. Therefore, on the basis of prior study, further analysis is needed to determine the mathematical characteristics of Richards. While fitting distribution data, the parameter B of Richards is ,0. As a result, the following new expression formula of Richards function in Table 3 can be derived:

Model Development
If While fitting diameter distributions [15], the empirical values of parameter B in Richards function are,23, then ln ({B).0, besides, parameter k is .0, it can be concluded that parameter q is .0, and parameter p also is .0. Because parameter m is .1, parameter r is ,0. Undoubtedly, Eq. (2) will still have high simulation precision, and it is obvious and important that the parameters of Eq. (2) may have the same theoretical meaning as three-parametric Weibull equation. Namely, parameter q may be the location parameter, parameter p may be the scale parameter, and parameter r may be the shape parameter. Note, whether the parameter q of Eq. (2) is.x or not, problems that the equation has no meaning or parameters have difficulty to converge will not appear. We call Eq.
(2) ''R distribution''. The probability density function of R distribution for a random variable x is Eqs. Accordingly, R distribution has a good application prospect in the field of diameter distribution. Our database, including 150 cumulative diameter distributions, was used to fit the R distribution and Weibull equation, then the theoretical meaning of parameters of R distribution was analyzed by discussing the correlativity between the parameters estimates and the whole stand characteristics.
The skewness 6 kurtosis can describle the shape and modeling properties of distribution function while using for forest stand [12]. The shape of R distribution is herein evaluated in terms of its skewness 6kurtosis from one side. Skewness and kurtosis values of frequency distributions of original data and R distribution are calculated for every 150 stands. The mathematical formulas are: where x i is the midpoint value of diameter class k, and d d is the average value of DBH of a stand, F i is the frequency of diameter class k, s is standard deviation of DBH.

Model Evaluation
Model evaluation was made through three steps: (1) the parameters of R distribution and three-parametric Weibull equstion were estimated by using nonlinear regression method (NRM) or maximum likelihood estimates method (MLEM) for each of the 150 unthinned stands, (2) prediction functions between parameters and stand characteristics were built by using parameter prediction method (PPM) and parameter recovery method (PRM) for all 150 unthinned stands, (3) the modelling and prediction properties were tested by adopting the 159 independent stands.
For R distribution, parameter estimates for p (p ) and c (c _ ) were obtained by employing NRM and MLEM. Because the location parameter a of the threeparametric Weibull function was often viewed as the minimum observed diameter or its multiple [32], parameter a was defined as the lower limit of the minimum diameter class when MLEM was adopted.
PPM and PRM were used to build stand-level diameter distribution models [5,[35][36][37][38]. PPM means direct regression of parameters from stand characteristics, and PRM means to recalculate parameters after estimating percentiles or some key point on equation curves from stand characteristics. Some parameters of R distribution and Weibull were regressed against stand characteristics which included stand age, site index, planting density, average height, dominant height and quadratic mean DBH (D g ), some parameters were derived from the moments of the diameter distribution, which were themselves estimated from stand diameter values in a stand, n is the total number of trees in the stand.
For R distribution, when PPM and PRM are simultaneously adopted to predict diameter distributions of stands used for model testing, based on preliminary graphical and correlation analyses, the parameter prediction equations for estimates for p _ , q _ and r _ are given by Eqs. (4), (5) and (6), respectively where x 1 ,x 2 ,x 3 ,x 4 ,x 5 ,x 6 are the stand characteristics, and I _ x is the estimate of abscissa of inflection point of stand diameter distribution, which can be derived by When only PRM is the adopted method, besides Eq. (6), Eqs. (8) and (9) can be used to recover the parameters where D  7).
For the three-parametric Weibull function, the PPM and PRM are simultaneously adopted to predict stand parameters. Parameters a _ and b _ can be estimated as parameters p _ and q _ of R distribution. Parameter c _ can be solved from the cumulative distribution function of Weibull by Eq. (10).
The reason that the 0.5 percentile was adopted is that this percentile was near the inflection point of the distribution curve [15]. The prediction of D To further test the modelling and prediction properties of R distribution and Weibull, another database comprised of 159 diameter frequency distributions was used to test the models.

Comparison of the Models
The application effect of R distribution and three-parametric Weibull function was examined by comparing the residual sum of square (RSS) and coefficient of determination (R 2 ). The residual sum of square and coefficient of determination were respectively calculated as where obs k and est k are the observed and predicted diameter frequency for diameter class k, and n is the number of diameter classes in a sample stand.
We also used the Kolmogorov-Smirnov test to test the goodness of fit of distribution functions. SAS 9.1 and EXCEL 2003 was adopted for parameter estimation and model evaluation.

Results and Discussion
Modelling Results of R Distribution and Three-Parametric Weibull Function Table 4 shows the mathematical characteristics of parameters and statistics for R distribution and the three-parametric Weibull function derived from the 150 plot measurements. For R distribution, the estimates of parameters p, q and r are .0, .0 and ,0, respectively, which is identical to the empirical distribution range mentioned in Equation (2). The floating range, mean and standard deviation of parameter p or q shows that these two parameters are both stable. The total RSS from all stands of R distribution (0.1913) was smaller than Weibull distribution (0.2069), and the mean of R 2 (R 2 = 0.9990) was slightly higher than Weibull distribution (R 2 = 0.9988) from every stand while NRM was applied. Although the precision of the two models are both high, R distribution presents a more accurate simulation than the three-parametric Weibull function (Table 4). For the threeparametric Weibull function, the precision of NRM is far higher than MLEM in the view of total RSS and R 2 ( Table 4).
The skewness 6 kurtosis of the original data and R distribution is depicted in Figure 1. Each circle dot represents a stand. The distribution range of skewness is similar between the original data and R distribution. Most stands have a negative skewness, which is different from inverted J distribution that often happens for the natural stand with positive skewness [2]. The kurtosis obtained from R distribution are a little smaller than the original data for some stands, but the value range of kurtosis is same as the actual values. The results mean the shapes modeled by R distribution are undistorted and multiple. Figure 2 shows both the actual data distribution (histogram), the estimated R distribution shapes (bold solid line), the estimated Weibull distribution obtained from NRM (solid line) and the estimated Weibull distribution obtained from MLEM (dashed line) for a selection of stands from different planting densities, stand ages and quadratic mean DBH. In general R distribution and Weibull distribution from NRM both have provided a good fit for all the stands analyzed, and Weibull distribution from MLEM provided a relatively bad fit. In comparison, R distribution is more stable than Weibull distribution from NRM, which can be seen from the application of the two distributions to stand 2, 4.

Theoretical Meaning of Parameters of R Distribution
The modelling accuracy and theoretical meanings of parameters are the two most important indexes used to judge whether a function is suitable to modelling diameter distributions and the two indexes are complementary. It is known that R distribution has high modelling precision (Table 4). However, do the parameters of R distribution have good theoretical interpretation? This question can be answered by discussing the relationship between parameters of R distribution and the basic stand characteristics such as stand age, stand density, site index and quadratic mean DBH. Figure 3 shows the relationship of parameter p of R distribution to stand age and quadratic mean DBH. Parameter p increased with increasing stand age and quadratic mean DBH. Liu et al. (2004) [38] also found that the relationship between scale parameter b of Weibull distribution and stand age was positive. The relationship of parameter p and stand age and quadratic mean DBH were well approximated by the brief second polynomial. The resultant parameter prediction equation for predicting p is given by Eq. (11).   The Theoretical Meaning of Parameters q of R Distribution

Theoretical Meaning of Parameters p of R Distribution
The relationship of parameter q of R distribution to stand age, stand density, site index and quadratic mean DBH is illustrated in Figure 4. Parameter q increased with increasing stand age, site index and quadratic mean DBH, while it decreased with increasing stand density. The coefficients of determination (R 2 ) of the second polynomial between parameter q and stand age, stand density, site index and quadratic mean DBH were respectively 0.2477, 0.5843, 0.3011 and 0.7886. It is obvious that quadratic mean DBH has the biggest effect on parameter q among these stand characteristics. The parameter prediction equation for predicting q is given by Eq. (12).  The coefficients of Eq. (12) were significant at the 0.0001 significance level. From Figure 4, it can be known that the positive or negative relativities of parameter q and the four stand characteristics can reasonably interpret the theoretical meaning of parameter q as a location parameter of diameter distribution. For Weibull distribution, the location parameter a also increased with increasing quadratic mean dbh and stand age [34,39]. However we should know the fact that the Weibull location parameter must be smaller than observed diameter in a stand. Based on our experience with the data set and other data sets, the convergence and precision of Weibull function are sensitive to the location parameter a. Some researchers regarded parameter a as minmum diameter in a stand or its times [4,340,41], such as 1/3, 1/2, and 1. And it is difficult to get a common sense. However, for R distribution, there is no strictly restrictive condition about location parameter q (eq. 2). According to the form of R distributin function, parameter q can be regarded as the location parameter of R distribution, which may make R distribution become a new promising diameter distribution.
In previous studies, Duan et al [42] discovered Richards function was suitable for modelling diameter distribution, but did not realize the underlying relationship between parameter B and k in Richards function. They thought parameter B had poor theoretical interpretation. For R distribution, parameter q, composed of parameters B and k, obviously has good theoretical meaning, and proved easy to converge. This might lead to the use of R distribution as a new diameter distribution.
Theoretical Meaning of Parameters r of R Distribution Figure 5 shows the relationship of parameter r of R distribution to stand age and quadratic mean DBH. Parameter r decreased with increasing stand age and quadratic mean DBH. The result of analysis of variance showed that parameter r had significant relativity with stand age and quadratic mean DBH at the 0.01 significance level. The coefficients of determination (R 2 ) of the second polynomial between parameter r and stand age and quadratic mean DBH were 0.0562 and 0.0706, respectively.
For a sigmoid curve, the location of inflection point decides its shape at some extent. Parameters p, q and r together decide the abscissa of inflection point of R distribution ( Table 3). As shown in where D g refers to quadratic mean DBH. The value of ((r{1)=r) r , being the ordinate of inflection point of R distribution, is decided by parameter r. Figure 7 shows the relationship of ordinate of inflection point of R distribution to stand age. The ordinate of inflection point decreased with increasing stand age. The coefficient of determination was 0.2402, which means that the ordinate of inflection point of R distribution can show the change of distribution shape with age in some extent. Besides, the potential relationship between the ordinate of inflection point of R distribution and its skewness is showed in Figure 8. It is found that the ordinate of inflection point of R distribution has significant relativity with its skewness and kurtosis, the linear coefficient of determination are 0.4483 and 0.2824 respectively. This phenomenon shows that the size of inflection point of R distribution can report the shape of stands in some extent.
Due to the act of parameter r in inflection point of R distribution, it can be concluded that parameter r has a deep relationship with the shape of R distribution. Besides, the significant correlation between abscissa and ordinate of inflection point of R distribution and stand characteristics indirectly shows that parameter r is flexible and theoretical. Accordingly, parameter r may be considered as the shape parameter of R distribution.

Distribution of Inflection Point of R Distribution
R distribution has a flexible inflection point. The variation range of the ordinates of inflection points is 0.3787,0.6436, and 91.33 percent of them distribute in the range of 0.4,0.6. Therefore, we can conclude that the main distribution range of inflection points for the cumulative diameter distribution of stands was 0.4,0.6.

Functions for the Estimation of Stand Parameters
Stepwise regression analysis was applied to build the relationship between the parameters of two distributions and the stand characteristics composed of stand age, planting density, site index and so on at a 0.5 risk level. The regressions are on the 150 estimation stands. When nonlinear regression method were adopted, the resultant parameter prediction equations for predicting p   Obviously, because of weak relativity, parameter r _ of R distribution was not adaptive to be directly predicted by stand characteristics, which could be indirectly obtained through Eqs. (6) or (13). For parameter a _ , b _ of Weibull distribution, the coefficients of determination of second degree polynomials between the other parameter and above-mentioned stand variable were all smaller than those of Eqs. (19), (20), (21), and (22), respectively. Therefore, the five equations were adopted to predict the unknown distribution parameters.
Note that for the three-parametric Weibull function, regardless of whether the nonlinear regression method or maximum likelihood estimates method is used, c _ always has significant relativity with planting density and stand age at the 0.0001 significance level. However, under the consideration of relativity, parameter c _ is obtained through Eqs. (21) and (22) or Eqs. (25), and (26). It can be found that these functions have different stand characteristics and exponentia, which are mainly concluded by the analysis of theoretical meaning of parameters and the needs of brief and high precision of models. In practice, while considering parameter prediction method, the parameter prediction models that have both acceptable high precision and theoretical meaning are firstly selected, which often include single stand characteristics, then for parameters with high relativity with some stand characteristics, stepwise regression analysis can be applied. While some or single parameter in distribution model with weak relativity with stand characteristics, parameter recovery method can be adopted. Of course, parameter recovery method can be applied in other conditions. Based on the above-mentioned parameter  Table 5. The regression expressions and coefficients of determination (R 2 ) between parameters or diameters at the key points and stand characteristics for R distribution and three-parametric Weibull. prediction equations that have a high coefficient of determination, parameters of distribution models can be evaluated by introducing the related stand characteristics into these equations.

Goodness-of-fit for Estimating Diameter Distribution
Data from 159 evaluation subplots provide an opportunity to analyze and compare the accuracy of the R distribution and threeparametric Weibull function. In these calculations, two fitting methods (nonlinear regression method and maximum likelihood estimates method) and two parameter estimation methods (PRM and the combination of PPM and PRM) were used.
It was encouraging that R distribution was found to have lower RSS and higher non-rejection rate than the three-parametric Weibull function (e.g. Table. 6). Figure 9 shows the distribution of the residual sum of square (RSS) against quadratic mean DBH, from which we can directly compare the prediction effects of the   Table 6. doi:10.1371/journal.pone.0062605.g009 five methods. Method A has the highest precision among the five methods, which shows R distribution can accurately estimate diameter distribution of most future stands using the combination of PPM and PRM under the condition that only two stand characteristics are known (e.g. Table. 6). For method A, the result of the Kolmogorov -Smirnov test showed the null hypothesis that the observed and fitted distributions that are the same cannot be rejected for 128 out of 159 stands (80.50%). However, in view of the practical application of distribution models to stand management, method B, based on R distribution, would be a better option because of its dependence on stand characteristics related to stand density and site quality and because its non-rejection rate reached 73.58% for the total test data. Methods B and D used the same parameter recovery method and numbers of stand characteristics and had a nearly equal non-rejection rate, which showed R distribution had a distribution function equally as good as threeparametric Weibull function. Furthermore, because the iterated function of R distribution was easier to converge than threeparametric Weibull function when using nonlinear regression method to predict parameters, R distribution would have a wide application prospect. For Chinese fir plantations, the non-rejection rate of unthinned stands and thinned stands did not have obvious differences when methods A, B, C and D were adopted to predict the diameter distributions of future stands (e.g. Table. 6).
Additionally, for the three-parametric Weibull function, the nonlinear regression method is a more effective approach than the maximum likelihood estimates method in the estimation system of stand diameter distribution (e.g. Table. 6).

Conclusion
Based on analysis of the disadvantages of the three-parametric Weibull function, this study develops a promising distribution function (R distribution), which is a new and essential exploration in the study of parametric methods. We conclude that: (1) R distribution has a more accurate simulation than the threeparametric Weibull function while modelling diameter distributions of Chinese fir plantations; (2) the parameters p, q and r of R distribution proved to be its scale, location and shape parameters, and have a deep relationship with stand characteristics; this means the parameters of R distribution have good theoretical interpretation; (3) the main distribution range of inflection points for the cumulative diameter distribution of Chinese fir plantations was 0.4,0.6; (4) the goodness-of-fit test showed the diameter distributions of unknown stands can be accurately estimated by applying R distribution and with regards to modelling precision and theoretical interpretation, method B (table 6) may be the most suitable choice due to its good convergence, high precision and including multiple stand characteristics.

Author Contributions
Conceived and designed the experiments: A-GD J-GZ. Performed the experiments: A-GD. Analyzed the data: A-GD. Contributed reagents/ materials/analysis tools: A-GD X-QZ C-YH. Wrote the paper: A-GD.