Modelling diameter distribution of Tetraclinis articulata in Tunisia using normal and Weibull distributions with parameters depending on stand variables

Abstract

The objective of this study was to evaluate the effectiveness of both Normal and two-parameter Weibull distributions in describing diameter distribution of Tetraclinis articulata stands in north-east Tunisia. The parameters of the Weibull function were estimated using the moments method and maximum likelihood approaches. The data used in this study came from temporary plots. The three diameter distribution models were compared firstly by estimating the parameters of the distribution directly from individual tree measurements taken in each plot (parameter estimation method), and secondly by predicting the same parameters from stand variables (parameter prediction method). The comparison was based on bias, mean absolute error, mean square error and the Reynolds’ index error (as a percentage). On the basis of the parameter estimation method, the Normal distribution gave slightly better results, whereas the Weibull distribution with the maximum likelihood approach gave the best results for the parameter prediction method. Hence, in the latter case, the Weibull distribution with the maximum likelihood approach appears to be the most suitable to estimate the parameters for reducing the different comparison criteria for the distribution of trees by diameter class in Tetraclinis articulata forests in Tunisia.

Introduction

Thuya (Tetraclinis articulata [Vahl.] Mast.) is a semiarid Mediterranean forest species which extends discontinuously from North Africa to south-western Europe. The natural character and the limited geographic distribution area of the species make it of particular ecological interest in terms of biodiversity conservation. Thuya forests provide highly valuable services such as soil erosion control, biodiversity conservation and CO₂ fixation, and are suitable for afforestation programs in arid or semiarid environments, where few species are able to grow ([48]), and areas with severely eroded soils ([19]).

In addition, these stands are of potential economic and social interest for the rural populations, not only in Tunisia where they cover 33000 ha ([16]), but also in Algeria and Morocco, covering more than 1000000 ha in Maghreb ([12]). Thuya timber is widely used in construction, handicraft and cabinetmaking due to its hardness, veneer, durability and fragrance. A secondary traditional use of this species is resin production to obtain sandarac gum. However, overexploitation for timber, irrational resin tapping, grazing pressure and uncontrolled fires have resulted in highly degraded forests ([12], [18]). The extinction of this species has only been avoided by its resprouting capacity. The current degraded state has left these forests devoid of large trees, so the main production presently comprises small pieces of wood for fuel and fencing ([17]). However, the profitability of these forests could be improved by exploiting other uses of the wood such as in industry, production of decorative objects, etc.

Improving the profitability of wood and resin production, along with the need to preserve and improve Tunisian thuya stands, which play a key ecological role in arid or semiarid environments, necessitates the development of models for quantitatively assessing the production of wood at different forest sites and under different management schedules to ensure a sustained yield. A first step towards this objective was the development of tree-level models relating common variables measured during forest inventories, such as breast height diameter, with relevant tree attributes like total tree height, crown diameter, height to crown base, stem taper curve and tree volume ([8]). The diameter-growth dynamics of the trees was studied by developing a distance-independent individual tree diameter increment model ([46]). A stand-level growth and yield model has recently been developed consisting of a site index sub-model and a system of compatible stand attribute equations presented in the form of a stand density management diagram ([47]). The next step would be to characterize the stand diameter structure.

Fonseca et al. ([21]) state that detailed information on the variability of tree diameters in natural and altered stands is key to sustainable forest management and is needed to assess carbon sequestration in forest ecosystems. Tree diameter distributions play an important role in stand management. For example, data regarding the number of trees per diameter class in a stand where a given silvicultural treatment has been applied are of particular interest because the diameter size will determine the industrial use of the wood and thus the price of the different products ([24]). Diameter distribution data at stand level also provides a more scientific basis for forest managers to choose which trees should be harvested ([45]).

Different types of parametric density functions have been used to describe tree diameter distribution in forest stands, including the Normal ([37]), Log-normal ([5], [37]), Gamma ([39]), Beta ([13], [42]), Johnson’s SB ([27], [42], [21]) and the Weibull distribution ([2], [42]). Mateus & Tomé ([33]) used Johnson’s distribution for modeling the diameter distribution of first-rotation eucalyptus plantations in Portugal. Gorgoso et al. ([25]) compared the accuracy of the Weibull, Johnson’s SB and beta distributions for modeling even-aged stands of Pinus pinaster Ait., Pinus radiata D. Don and Pinus sylvestris L. in northwest Spain, and Binoti et al. ([4]) evaluated the effectiveness of fatigue life, Frechet, Gamma, Generalized Gamma, Generalized Logistic, Log-logistic, Nakagami, Beta, Burr, Dagum, Weibull and Hyperbolic distributions for describing diameter distribution Tectona grandis L. stands subjected to thinning treatments at different ages in Brazil. Sghaier & Palm ([45]) used Pearson’s type I function for modeling diameter distribution of Pinus halepensis Mill. stands in Tunisia.

The objective of this study is to model the diameter distribution of the main Tetraclinis articulata stands in Tunisia. To achieve this objective, we compared the Normal distribution, given its ease of use, and the two-parameter Weibull distribution, since it has been shown to be the most suitable in many cases ([7], [34]).

Material and ethods

Study area and data collected

The natural Tetraclinis articulata stands at Jbel Lattrech forest in north-eastern Tunisia cover an area of approximately 4.000 hectares and are essentially mono specific and homogeneous. The climate is semi-arid, with an annual average precipitation of 393 mm, of which 95 % of rain falls between September and May. The average minimum temperatures for the coldest month (January) and maximum for the hottest month (August) are 8.4 and 30.5 °C, respectively ([3]). The study zone is located on Oligocene sandstones where there are groupings of Tetraclinis articulata - Lavandula stoechas, dominated by Cistus monspelliensis L., Genista aspalathoides Lam., Erica arborea L., Brachypodium ramosum Roemer et Schultes, Avena bromoides Gouan and Lupinis angustifolius L. ([36]).

The data used to develop the diameter distribution model were collected in 2009 from 50 temporary circular sample plots of variable radius and established such that they cover a wide range of age, density and site quality combinations. The size of the plots ranged from 88 to 835 m², depending on the stand density, so as to include at least 40 measured trees per plot. The diameter at breast height of all trees in each plot (down to a minimum diameter of 5 cm) was measured to the nearest 0.1 cm. Additionally, the heights of a random sample of 11 trees per plot (including the 4 largest trees) were measured to estimate mean height. Dominant height was calculated from the four largest trees per plot, and age was estimated from stem analysis on fallen trees by counting the number of rings at 0.30 m (see [8] for more details on data collection). Tab. 1 shows the main stand variables measured in the 50 plots.

Tab. 1 - Summary of main stand variables for the sample used for modelling ([47]).

Variables	Mean	STD	Min	Max
Density (stems ha-1)	1860	865	479	4533
Age (years)	57	10	36	77
Site index (m)	5.06	1.08	2.83	7.78
Basal area (m2 ha-1)	12.90	6.39	2.32	27.74
Quadratic mean diameter (cm)	9.41	1.55	6.32	12.95
Dominant diameter (cm)	12.76	2.32	7.85	17.43
Dominant height (m)	6.52	1.38	4.20	11.00
Mean height (m)	5.46	1.09	3.61	9.17
Volume (m3 ha-1)	36.14	22.13	5.16	111.61

Diameter distribution unctions

We tested two distribution functions: the Normal function, known for its ease of use ([29]), and the Weibull function, most frequently used because of its flexibility, simplicity and existence of a closed form solution for its CDF ([10], [35], [41], [24]).

Normal distribution

The equation for the normal probability density function is (eqn. 1):

(1)f(x)=σπexp[−(x−mσ)]" role="presentation"> f(x)=σπ−−√exp[−(x−mσ)](1) (1)f(x)=1σ2πexp[−12(x−mσ)2]

where x is the random variable, and m and σ are its arithmetic mean and standard deviation, respectively.

The cumulative probability distribution of the normal function is (eqn. 2):

(2)F(x)=∫∞xf(x)dx" role="presentation"> F(x)=∫x∞f(x)dx(2) (2)F(x)=∫∞xf(x)dx

The two parameters which define the normal distribution (arithmetic mean and the standard deviation) can be estimated by using the following expressions ([15] - eqn. 3, eqn. 4):

(3)m^=n∑i=nxi" role="presentation"> m^=n∑i=nxi(3) (3)m^=1n∑i=1nxi

(4)σ^=n−∑i=n(xi−m^)" role="presentation"> σ^=n−∑i=n(xi−m^)−−−−−−−−−−−−−−−−√(4) (4)σ^=1n−1∑i=1n(xi−m^)2

where n indicates the number of trees per plot and x_i (cm) the diameter at breast height of each tree.

Weibull distribution

The variety of shapes that can be taken by the Weibull distribution makes it particularly suitable for representing tree diameter distribution. The equation for the Weibull probability density function is as follows (eqn. 5):

(5)f(x)=cb(x−ab)c−exp[−(x−ab)c]" role="presentation"> f(x)=cb(x−ab)c−exp[−(x−ab)c](5) (5)f(x)=cb(x−ab)c−1exp[−(x−ab)c]

for x ≥ a, a ≥ > > or if not. The cumulative probability distribution of the Weibull function is (eqn. 6):

(6)F(x)=−exp[−(x−ab)c]" role="presentation"> F(x)=−exp[−(x−ab)c](6) (6)F(x)=1−exp[−(x−ab)c]

for x > a. This function is defined by three parameters: a is the location parameter, b is the scale parameter and c is the shape parameter.

Parameter a, commonly termed the location parameter, identifies the lower bound of the diameter distribution. For fixed values of b and c, changes in the parameter a simply shift the entire distribution along the x-axis. Parameter b is the scale parameter, and point x = a + b corresponds approximately to the 63rd percentile of the distribution. Parameter c defines the shape of the Weibull distribution. When cc = 1, the Weibull distribution is reduced to the exponential distribution. When c > 1, the distribution is mound-shaped, approximately equivalent to the Normal distribution for c = 3.6. When c is between 1 and 3.6 the Weibull distribution is positively skewed, and is negatively skewed when c is greater than 3.6 ([2]).

The most frequently used methods to estimate the parameters of the Weibull distribution, are the maximum likelihood method ([28]) and the method of moments ([6]). Parameter a can be considered the smallest possible diameter in the stand and thus should be between zero and the minimum observed value in some cases ([2]). In some studies, the parameter a is arbitrarily fixed at 0.5 d_min ([30]) or at zero ([24], [44]), thus reducing the function to a two-parameter Weibull distribution which is easier to model and provides similar results to those of the three-parameter Weibull function, at least for even-aged, single-species stands ([34], [32]).

The two-parameter Weibull function (with the location parameter a = 0) was considered in this study as follows (eqn. 7):

(7)f(x)=cb(xb)c−exp[−(xb)c]" role="presentation"> f(x)=cb(xb)c−exp[−(xb)c](7) (7)f(x)=cb(xb)c−1exp[−(xb)c]

where x is the random variable, b the scale parameter and c the shape parameter that controls the skewness.

Maximum likelihood estimator (MLE)

The maximum likelihood method is a commonly used procedure for the Weibull distribution in forestry because it has certain desirable properties ([30]). Estimation of the parameters using maximum likelihood has been found to produce consistently better goodness-of-fit statistics compared to other methods, but it also puts the greatest demands on the computational resources ([11]). If we consider the Weibull probability density function given in eqn. 7, then the likelihood function (L) will be (eqn. 8):

(8)L(x,...,xn;b,c)=∏i=ncb(xib)c−exp[−(xib)c]" role="presentation"> L(x,...,xn;b,c)=∏i=ncb(xib)c−exp[−(xib)c](8) (8)L(x1,...,xn;b,c)=∏i=1ncb(xib)c−1exp[−(xib)c]

Taking the logarithms from eqn. 8, differentiating with respect to b and c respectively, and satisfying the following equations ([14], [38], [20] - eqn. 9, eqn. 10):

(9)c=[∑i=nxiclog(xi)∑i=nxic−n∑i=nlog(xi)]−" role="presentation"> c=[∑ni=xcilog(xi)∑ni=xci−n∑i=nlog(xi)]−(9) (9)c=[∑i=1nxiclog(xi)∑i=1nxic−1n∑i=1nlog(xi)]−1

(10)b=(n∑i=nxic)c−" role="presentation"> b=(n∑i=nxci)c−(10) (10)b=(1n∑i=1nxic)c−1

where n equals the number of sample observations in a Weibull distribution and x_i the diameter of each tree. The value of c must be obtained from eqn. 9 by using standard iterative procedures and then it is used in eqn. 10 to obtain b.

Method of moments (MOM)

The method of moments is another technique commonly used for parameter estimation. In the Weibull distribution, the k moment readily follows from eqn. 7 and is ([30] - eqn. 11):

(11)mk=(b)k/cΓ(+kc)" role="presentation"> mk=(b)k/cΓ(+kc)(11) (11)mk=(1b)k/cΓ(1+kc)

where Γ (Gamma function) is Γ(s) = ∫ x^s-1 e^-x dx (s > 0). Then from eqn. 11, we can find the first and the second moment as follows (eqn. 12, eqn. 13):

(12)m=μ=(b)/cΓ(+c)" role="presentation"> m=μ=(b)/cΓ(+c)(12) (12)m1=μ=(1b)1/cΓ(1+1c)

(13)m=μ+σ=(b)/cΓ(+c)" role="presentation"> m=μ+σ=(b)/cΓ(+c)(13) (13)m2=μ2+σ2=(1b)2/cΓ(1+2c)

which gives (eqn. 14):

(14)σ=m−μ=(b)/c[Γ(+c)−Γ(+c)]" role="presentation"> σ=m−μ=(b)/c[Γ(+c)−Γ(+c)](14) (14)σ2=m2−μ2=(1b)2/c[Γ(1+2c)−Γ2(1+1c)]

where σ² is the variance of tree diameters in a plot, and m₁, m₂ are the arithmetic and quadratic mean diameter in a plot, respectively.

When σ² is divided by the square of m₁, the expression for obtaining c is (eqn. 15):

(15)σμ=Γ(+c)−Γ(+c)Γ(+c)=σ=d¯Γ(+c)⋅[Γ(+c)−Γ(+c)]" role="presentation"> σμ=Γ(+c)−Γ(+c)Γ(+c)=σ=d¯Γ(+c)⋅[Γ(+c)−Γ(+c)](15) (15)σ2μ2=Γ(1+2c)−Γ2(1+1c)Γ2(1+1c)=σ2=d¯2Γ2(1+1c)⋅[Γ(1+2c)−Γ2(1+1c)]

In order to estimate b and c, we need to calculate the arithmetic mean diameter d and the variance σ² of the observed distribution and obtain the estimator of c in eqn. 15. eqn. 15 was resolved by an iterative procedure. When the value of the location parameter (a) is zero, the scale parameter (b) can then be calculated directly using the following equation ([24], [44] - eqn. 16):

(16)b=d¯Γ(+c)" role="presentation"> b=d¯Γ(+c)(16) (16)b=d¯Γ(1+1c)

where d is the arithmetic mean diameter.

Method comparison

To determine whether the adjusted theoretical distribution provides a good representation of what is observed, we used the chi-square test (χ²). This test requires some extreme classes to be grouped when there are an insufficient number of observations ([15]). For this reason, the following goodness-of-fit statistics were computed for each method: the mean bias, which reflects the deviation of the model with respect to observed values and the root mean-square error (RMSE), which is a measure of the precision of the estimates.

For the calculation of these statistics, the observed value Y_i and the theoretical value Y_i were calculated for each diameter class per plot (2 cm diameter class). In addition, the error index ([43]) was calculated (eqn. 17):

(17)EIi=∑i=k|Yik−Y^ik|" role="presentation"> EIi=∑i=k∣∣Yik−Y^ik∣∣(17) (17)EIi=∑i=1k|Yik−Y^ik|

where Y_i and Y_i are the observed and predicted number of trees in diameter class k for the i-th plot. The sum includes all diameter classes in the i-th plot. In this error index, we can multiply the absolute difference between the observed and predicted number of trees per diameter class by the volume, using an equation developed by Sghaier et al. ([47]): v = 2.182 · 10^-2 · d^2.89, R² = 0.917. The diameter used in this equation corresponds to the middle diameter class. In this case (eqn. 18):

(18)EIi=∑i=k(|Yik−Y^ik|⋅vk)" role="presentation"> EIi=∑i=k(∣∣Yik−Y^ik∣∣⋅vk)(18) (18)EIi=∑i=1k(|Yik−Y^ik|⋅vk)

This index can also be expressed relatively as a percentage by dividing it by the total volume of the observed distribution in the plot ([29] - eqn. 19).

(19)EIi′=EIi∑i=kYivi⋅" role="presentation"> EI′i=EIi∑ki=Yivi⋅(19) (19)EIi′=EIi∑i=1kYivi⋅100

The latter error index was systematically calculated for each plot in this study.

Parameter prediction: distribution parameters with stand variables

In order to predict the parameters of the distribution functions directly from the stand variables using simple regression models, correlation analyses were carried out by estimating Pearson’s correlation coefficients between the estimated parameters corresponding to each plot and different stand variables. The stand variables chosen were: age (A), density (N), dominant height (H_d), mean height (H_m), quadratic mean diameter (d_g), arithmetic mean diameter (d), 25% percentile (P₂₅), median (P₅₀), 75% percentile (P₇₅) of the stand, and natural logarithm transformations for the quotient of the 25% percentile, the 50% percentile, the 75% percentile and the mean diameter by the quadratic mean diameter (LP₂₅, LP₅₀, LP₇₅, and Ld, respectively).

To predict the fitting parameters from the yield table using the Normal and Weibull distribution, the mean diameter and the variance must be known (first and second order moments of the distribution, respectively). The variance can be directly obtained from the arithmetic and the quadratic mean diameter ([22]) by the following expression (eqn. 20):

(20)σ=dg−d¯" role="presentation"> σ=dg−d¯(20) (20)σ2=dg2−d¯2

The mean diameter can be modeled with eqn. 21, which ensures the prediction of dÌis lower than the quadratic mean diameter (d_g - eqn. 21):

(21)d¯=dg−eXβ" role="presentation"> d¯=dg−eXβ(21) (21)d¯=dg−eXβ

where X is a vector of independent stand variables in a fixed instant and β is the vector of parameters to be estimated ([24], [44]).

To select the model, the bias, the mean square error and the adjusted determination coefficient were computed.

Results

Tab. 2 shows summary descriptive statistics for the estimated parameters.

Tab. 2 - Summary statistics of the estimated parameters for the studied distributions (50 plots). (Sigma): σ; (predicted mean): μ.

Distribution	Method	Parameter	Mean	Minimum	Maximum	CV%
Weibull	MLE	b	10.01	6.680	13.794	16.5
		c	4.134	3.137	7.087	18.6
	MOM	b	9.986	6.657	13.777	16.6
		c	4.418	3.077	7.410	20.8
Normal	-	Predicted mean	9.087	6.240	12.550	16.1
		Sigma	2.427	1.011	3.862	26.3

In order to adjust the models to predict the parameters of the theoretical distributions using the stand variables, Pearson’s linear correlations for estimated parameters were calculated (Tab. 3).

Tab. 3 - Correlations between stand variables and parameters estimated with the three studied distributions. (***): P < 0.001; (**): P < 0.01; (*) P < 0.05; (A): age (years); N: density (trees ha^-1); H_d: dominant height (m); H_m: mean height (m); d_g: quadratic mean diameter (cm); (dÌ ): mean diameter (cm); P₂₅: 25% percentile (cm); P₅₀: 50% percentile (cm); P₇₅: 75% percentile (cm); LP₂₅: ln(P₂₅/d_g); LP₅₀: ln(P₅₀/d_g); LP₇₅: ln(P₇₅/d_g); Ld=ln(d/d_g); ln: natural logarithm.

Variable	Normal distribution	Weibull distribution
		Maximum Likelihood		Moments
		b	c	b	c
A	0.088	0.218	0.176	0.221	0.110
N	-0.054	-0.170	-0.066	-0.174	-0.081
H d	0.665***	0.667***	-0.368**	0.670***	-0.464***
H m	0.720***	0.762***	-0.335*	0.763***	-0.456***
d g	0.872***	0.999***	-0.380**	0.999***	-0.511***
d	0.844***	0.998***	-0.336*	0.998***	-0.469**
P 25	0.648***	0.923***	-0.109	0.923***	-0.226
P 50	0.768***	0.970***	-0.231	0.971***	-0.367**
P 75	0.857***	0.976***	-0.359*	0.976***	-0.505***
LP 25	-0.617***	-0.253	0.742***	-0.254	0.783***
LP 50	-0.370**	-0.064	0.602***	-0.061	0.582***
LP 75	0.203	0.203	-0.042	0.206	-0.155
Ld	-0.811***	-0.463***	0.884***	-0.459***	0.908***

The highest Pearson correlation between the parameters of Normal and Weibull distribution and the stand variables was obtained with d_g and its natural logarithm transformations. The values obtained for Pearson’s linear correlation between the Ld parameter and the two estimated parameters of the Weibull distribution using the Maximum Likelihood and moments approaches were -0.463 and -0.459 for parameter b and 0.884 and 0.908 for parameter c. In the case of the shape parameter (c), the highest values obtained for Pearson’s linear correlation were for the variable Ld (Tab. 3). The values obtained for the coefficients of determination were 91% and 96% for models that relate the parameter c to the variable Ld (Tab. S1 in Appendix 1).

As regards the variance of the Normal distribution and the scale parameter (b) for the Weibull distribution, the results also show good linear relationships between these parameters and the measured stand variables. It appears that there is a logarithmic relation between the variance and quadratic mean diameter (d_g). The coefficient of determination was 78% (Tab. S1 in Appendix 1).

Note that the mean diameter is needed to estimate the shape parameter (c) in both parameter prediction methods of the Weibull distribution as well as in the Normal distribution. Since the yield tables generally only give the quadratic mean diameter (d_g), an equation was employed to predict the mean diameter from the quadratic mean diameter (d_g), the mean height (H_m) and the age (A) of the stand (Tab. S2 in Appendix 1).

According to the fitted prediction equation of the mean diameter (Tab. S2 in Appendix 1), the mean quadratic diameter (d_g), the mean height (H_m) and the age of the stand (A) must be known. Of these three parameters, only the mean height (H_m) cannot be directly obtained from yield tables and was therefore estimated by the fitting equation (eqn. 22) using the dominant height (H_d) disaggregated from yield tables as an independent variable (eqn. 22):

(22)Hm=⋅Hd" role="presentation"> Hm=⋅Hd(22) (22)Hm=0.8307⋅Hd

with R² = 0.989 and CV% = 11.03. The expected theoretical values for each plot are calculated using the three studied distributions with the two parameter prediction methods (estimated and prediction parameters). The observed, estimated and predicted values were grouped for each plot in 2 cm diameter classes. The results of the different comparison criteria are shown in Tab. 4.

Tab. 4 - Bias, mean absolute error, residual mean square error and mean error index in number of trees for the three compared distributions (Weibull MLE, Weibull MOM and Normal) and the two methods of construction (parameter estimation and parameter prediction).

Methods	Comparison criteria	Weibull		Normal
		Maximum Likelihood	Moments
Parameters estimation	Bias	0.0097	0.0068	0.0185
	MAE	1.7140	1.7145	1.6197
	RMSE	5.2173	5.3054	4.7749
	EI’ %	32.61	33.40	31.42
Parameters prediction	Bias	0.0164	0.0063	0.0220
	MAE	1.7180	1.7592	1.8958
	RMSE	5.3606	5.6220	6.1839
	EI’ %	32.25	34.03	35.13

In the parameter estimation, the Normal distribution is the most suitable for all the comparison criteria evaluated, with the exception of bias. Furthermore, the two parameter estimation approaches of the Weibull function gave similar results, although the maximum likelihood approach (MLE) was found to be slightly better.

In contrast to the parameter prediction in the first model construction approach (parameter estimation), the results obtained for parameter prediction revealed that the Weibull distribution was better for all the comparison criteria used. In this case, the maximum likelihood method also provided the best results, with the exception of the bias, for which the method of moments gave the lowest value.

In addition to comparing the distributions according to the Reynolds’ index error (EI %), which is calculated for each plot, an analysis of variance (ANOVA) was performed for each model construction method. The results obtained reveal that there is no significant difference between the three distributions in each of the two model construction methods (α = 0.75 for parameter estimation and α = 0.57 for parameter prediction, respectively).

The Weibull distribution, using the method of moments, occupied the same intermediate position with the exception of the bias, for which it gave the lowest value in both cases (estimated and prediction parameters). To analyze these results in greater depth, some additional graphs were produced. The first two, (Fig. 1 and Fig. 2), present the values of Bias and RMSE in each diameter class. Fig. 3 shows the mean value of the error index (EI %) in each density class.

Fig. 1 - Values of bias in number of trees in each diameter class obtained using two fitting methods of the two-parameter Weibull distribution and Normal distribution. (a) Parameter estimation; (b): parameter prediction.

Fig. 2 - Values of mean square error (RMSE) in number of trees for each diameter class obtained using the two fitting methods of the two-parameter Weibull distribution and Normal distribution. (a) Parameter estimation; (b) parameter prediction.

Fig. 3 - Mean value of error index EI % in each density class (trees ha-1). (a) Parameter estimation; (b) parameter prediction.

According to the values for bias (Fig. 1), all the distributions underestimated the frequency of trees in the smallest and largest diameter classes, and overestimated them for the 10 cm diameter class under both parameter construction methods. With the exception of the 8 cm diameter class, for which the Weibull maximum likelihood method gives the highest value, the distributions compared present a similar level of precision, though the Normal distribution is slightly better, especially for parameter prediction.

As regards the mean square error (Fig. 2), although there are slight differences in the smallest diameter classes (lower or equal to 10 cm), the distributions gave similar results for diameter classes of 12 cm and above, being the Weibull distributions slightly superior in the case of parameter prediction.

For the Reynolds error index (EI %), Fig. 3 shows the variation of the mean value of this criterion according to stand density class (trees ha^-1). For both model construction approaches, the Weibull distribution, using the maximum likelihood method to estimate the parameters, appears to give the most stable results of the three studied methods. In addition to its stability across all the density classes, this method also provides the most accurate results for the parameter prediction approach (Fig. 3).

According to the distribution of trees per diameter class in each of the 50 measured plots, three types of distribution shape were observed (symmetric distribution, right dissymmetric distribution and reverse-J distribution). In order to better appreciate the capacity of the adjusted models to predict the number of trees per diameter class at plot level, one representative plot per distribution type was selected. Graphs for estimated vs. observed values were constructed for a symmetric distribution (Fig. S1 in Appendix 1), a right skewed distribution (Fig. S2 in Appendix 1) and a reverse-J shaped distribution (Fig. S3 in Appendix 1) for the parameter estimation and parameter prediction methods. The three figures reveal that all the studied models behave similarly, though the Weibull function appears to be marginally better than the others, in particular for the parameter prediction approach.

The use of the model for predicting the distribution of trees by diameter class

The diameter distribution models can be used independently by measuring a number of stand variables, or together with the yield table, which provides information that allows to predict the parameters of the distribution employed. The example which is provided here uses information from the yield table elaborated as part of this study for modeling the growth and production associated with the management and conservation of Tetraclinis articulata forests in Tunisia ([47]). Tab. 5 presents the predicted Weibull PDF (probability density function) parameters for different ages of the best site quality (H_d = 7 m at 50 years of age) in the case of low stand density (Reineke stand density index: SDI

Tab. 5 - Values of Weibull PDF parameters at 20, 40, 60, 80 and 100 years of age for stands with SDI

Stand parameters					Weibull PDF parameters
Age	N	H m	D g	D m	b	c
20	1909	3.43	7.21	6.92	7.51	3.65
40	1400	5.81	10.1	9.70	10.84	3.65
60	1166	7.74	12.12	11.61	12.64	3.58
80	1024	9.33	13.67	13.07	13.83	3.53
100	925	10.26	14.91	14.31	14.69	3.64

Fig. 4 shows the curves associated with these distributions and Tab. 6 represents the distribution of trees (in number and as a percentage) per diameter class for various stand ages.

Fig. 4 - Curves of diametric distribution at 20, 40, 60, 80 and 100 years by diameter (DBH, cm) for stands with SDI [ less than ] 400 and SI =7 m (density of probability).

Discussion

Our Pearson’s correlations analysis results are very similar to the findings of Gorgoso et al. ([24]) for Betula alba stands in northwest Spain, using the two-parameter Weibull function, except for the parameter Ld, which expresses the logarithm of the mean diameter divided by the quadratic mean diameter. In addition, the good linear correlations which are obtained between the scale parameter (b) of the Weibull function and the stand characteristics are analogous to those observed by Campos & Leite ([9]), Leite et al. ([31]) and Binoti et al. ([4]). Linear models that relate b to the logarithm of the quadratic mean diameter (d_g) show a very high determination coefficient. Similar linear models for this relationship, with values of the adjusted determination coefficient close to 99%, were obtained by Alvarez González ([1]) for Pinus pinaster in Galicia and García-Guëmes et al. ([26]) for Pinus pinea stands in Valladolid (Spain).

The results obtained in this study reveal that the Normal distribution, though less flexible, provides results which are as satisfactory as those of the Weibull distribution. Lejeune ([29]) reported similar findings for Picea abies in Belgium. However, the Weibull distribution using the maximum likelihood method to estimate parameters appears to provide a slightly superior method in terms of RMSE, MAE and EI%, especially with the parameter prediction approach.

A number of authors have reported that the two-parameter Weibull function provides the most simple and accurate approach for modelling diameter distributions, including Gorgoso ([23]), and [24]) for birch stands (Betula alba L.) in north-western Spain, Maltamo et al. ([34]) for Pinus sylvestris and Picea abies (L.) Karst. stands in Finland, Alvarez González ([1]) for Pinus pinaster stands, and Condés ([14]) for different species in Spain. Furthermore, the two-parameter Weibull function has been used in a large number of related studies due to its flexibility and high correlation of its parameters with stand characteristics ([40], [4]). Therefore, given its flexibility and the very strong correlation between the PDF parameters and stand characteristics, the two-parameter Weibull distribution, using the maximum likelihood method, appears to be the most suitable for describing diameter distributions in this study.

Gorgoso et al. ([25]) and Zhang et al. ([49]) in their studies in north west Spain and the north-eastern part of North America, respectively, obtained better results with the 3-parameter Weibull distribution, using the maximum likelihood method rather than the moments approach.

Conclusion

In this study, a model has been developed for predicting the distribution of trees by diameter class using stand variables from Tetraclinis articulata forests in north-eastern Tunisia. The developed model, used in conjunction with the previously developed models for tree level attributes ([8]), the distance-independent individual tree diameter-increment model ([46]) and the stand based growth and yield model presented in Sghaier et al. ([47]), constitutes a useful tool for the sustainable management of Tetraclinis articulata forests in north-east Tunisia.

Acknowledgements

This study (Project Number: A/017275/08) was carried out within the framework of the bilateral co-operation between Tunisia (INRGREF) and Spain (INIA), and was financed by the Ministry for Scientific Research and Technology of Tunisia (MSRT) and the Spanish Agency of International Cooperation and Development (AECID), for which they are greatly acknowledged. We also thank Adam Collins for the English review of the manuscript.

References

(1) Alvarez González JG (1997). Estudio de distribuciones diamétricas en masas regulares de Pinus pinaster Ait. en Galicia. [Diameter distribution studies of Pinus pinaster Ait. in Galicia]. PhD thesis, Universidad Politécnica de Madrid, ETS de Ingenieros de Montes, Madrid, Spain, pp. 270. Gscholar

(2) Bailey RL, Dell TR (1973). Quantifying diameter distributions with the Weibull function. Forest Science 19: 97-104. Gscholar

(3) Ben Mansoura A, Garchi S (2001). Caractérisation de la croissance et de la régénération du Thuya par une technique modifiée de mesure de distances. [Characterization of the growth and the regeneration of the Thuya by a modified technique of distances measurement]. Les Annales de l’INRGREF, Numéro spécial 2001: 54-76. Gscholar

(4) Binoti DHB, Binoti MLMS, Leite HG, Fardin L, Oliveira JC (2012). Probability density functions for description of diameter distribution in thinned stands of Tectona grandis. Cerne, Lavras 18 (2): 185-196.

(5) Bliss CI, Reinker KA (1964). A log normal approach to diameter distributions in even-aged stands. Forest Science 10: 350-360. Gscholar

(6) Burk TE, Newberry JD (1984). A simple algorithm for moment-based recovery of Weibull distribution parameters. Forest Science 30 (2): 329-332. Online | Gscholar

(7) Borders BE, Souter RA, Bailey RI, Ware KD (1987). Percentile-based distributions characterize forest stands tables. Forest Science 33 (2): 570-576. Online | Gscholar

(8) Calama R, Sánchez-González M, Garchi S, Ammari Y, Cañellas I, Sghaier T (2012). Towards the sustainable management of thuya (Tetraclinis articulata (Vahl.) Mast.) forests in Tunisia: models for main tree attributes. Forest Systems 21 (2): 210-217.

(9) Campos JCC, Leite HG (2009). Mensuração florestal: perguntas e respostas. [Forest measurements: questions and answers] (3rd edn). Universidad Federal de Viçosa, Viçosa, MG, Brazil, pp. 548. Gscholar

(10) Cao QV (2004). Predicting parameters of a Weibull function for modelling diameter distributions. Forest Science 50: 682-685. Online | Gscholar

(11) Cao QV, McCarty SM (2006). News methods for estimating parameters of Weibull functions to characterize future diameter distributions in forest stands. In: Proceedings of the “13th Biennial Southern Silvicultural Research Conference”. General Technical Report SRS-92, Southern Research Station, USDA Forest Service, Asheville, NC, USA, pp. 338-340. Gscholar

(12) Charco J (1999). El bosque mediterráneo en el norte de África: biodiversidad y lucha contra la desertificación. [The Mediterranean forest in North Africa: biodiversity and fight against the desertification]. Ed. Mundo Árabe e Islam - AECI, Madrid, Spain, pp. 370. Gscholar

(13) Clutter JL, Bennett FA (1965). Diameter distribution in old-field slash pine plantations. Report no. 13, Georgia Forest Research Council, Macon, GA, USA, pp. 9. Gscholar

(14) Condés S (1997). Simulación de parcelas arboladas con datos del II Inventario Forestal Nacional [Simulation of plots hoisted with data of the II National Forest Inventory]. PhD thesis, Universidad Politécnica de Madrid, ETS de Ingenieros de Montes, Madrid, Spain, pp. 616. Gscholar

(15) Dagnelie P (1998). Statistique théorique et appliquée. Tome 1: Statistique descriptive et base de l’inférence statistique. [Descriptive statistics and basis of the statistical inference]. De Boeck, Bruxelles, Belgium, pp. 508. Gscholar

(16) DGF (1995). Résultats du premier inventaire forestier national en Tunisie [Results of the first national forest inventory in Tunisia]. Direction Générale des Forêts, Tunisie, pp. 88. Gscholar

(17) DREF (2002). Thuya: importance écologique et économique. [Thuya: Ecological and economic importance]. Terre et Vie 52: 4. Gscholar

(18) El-Mouridi M, Laurent T, Famiri A, Kabouchi B, Alméras T, Calchéra G, El Abid A, Ziani M, Gril J, Hakam A (2011). Caractérisation physique du bois de la loupe de Thuya (Tetraclinis articulata (Vahl) Masters). [Physical characterization of the magnifying glass wood of the Thuya (Tetraclinis articulata (Vahl) Masters)]. Physical Chemical News 59: 57-64. Online | Gscholar

(19) Esteve-Selma J, Martinez-Fernandez J, Hernández I, Montávez JP, Lopez JJ, Calvo JF, Robledano F (2010). Effects of climatic change on the distribution and conservation of Mediterranean forests: the case of Tetraclinis articulata in the Iberian Peninsula. Biodiversity and Conservation 19: 3809-3825.

(20) Eerikäinen K, Maltamo M (2003). A percentile based basal area diameter distribution model for predicting the stand development of Pinus kesiya plantations in Zambia and Zimbabwe. Forest Ecology and Management 172: 109-124.

(21) Fonseca T, Marques C, Parresol BR (2009). Describing maritime pine diameter distributions with Johnson’s SB distribution using a new all-parameter recovery approach. Forest Science 55 (4): 367-373. Online | Gscholar

(22) Frazier JR (1981). Compatible whole-stand and diameter distribution models for loblolly pine. PhD thesis, School of Forestry and Wildlife, Virginia Polytechnic Institute and State University, Blackburg, VA, USA, pp. 125. Gscholar

(23) Gorgoso JJ (2003). Caracterización de las distributions diametricas de Betula alba L. en Galicia. [Characterization of diameter distribution of Betula alba L. in Galicia]. PhD thesis, Universidad Santiago de Compostela, Chile, pp. 176. Gscholar

(24) Gorgoso JJ, Alvarez González JG, Rojo A, Grandas-Arias JA (2007). Modelling diameter distributions of Betula alba L. stands in northwest Spain with the two-parameter Weibull function. Investigation Agraria Sistemas Recursos Forestales 16 (2): 113-123.

(25) Gorgoso JJ, Rojo A, Camara-Obregon A, Dieguez-Aranda U (2012). A comparison of estimation methods for fitting Weibull, Johnson’s SB and beta functions to Pinus pinaster, Pinus radiata and Pinus sylvestris stands in northwest Spain. Forest Systems 21 (3): 446-459.

(26) García-Guëmes C, Cañadas N, Montero G (2002). Modelización de la distribución diamétrica de las masas de Pinus pinea L. de Valladolid (España) mediante la función de Weibull. [Modelling diameter distribution of Pinus pinea L. stands in Valladolid (Spain) by using Weibull function]. Investigation Agraria Sistemas Recursos Forestales 11: 264-282. Gscholar

(27) Hafley WL, Schreuder HT (1977). Statistical distributions for fitting diameter and height data in even-aged stands. Canadian Journal of Forest Research 4: 481-487.

(28) Jonhson NL, Kotz S (1970). Distribution in statistics: continuous univariate distributions. Houghton Mifflin Co, Boston, MA, USA, vol 2, pp. 306. Gscholar

(29) Lejeune P (1994). Construction d’un modèle de répartition des arbres par classes de grosseur pour les plantations d’épicéa commun (Picea abies L Karst) en Ardenne belge. [Construction of a distribution model of trees by diameter classes for the common spruce (Picea abies L Karst) plantations in Belgian Ardenne]. Annals of Forest Science 51: 53-65.

(30) Lei Y (2008). Evaluation of three methods for estimating the Weibull distribution parameters of Chinese pine (Pinus tabulaeformis). Journal of Forest Science 54 (12): 566-571. Online | Gscholar

(31) Leite HG, Binoti DHB, Guimarães DP, Silva MLM, Garcia SLR (2010). Avaliação do ajuste das funções Weibull e hiperbólica a dados de povoamentos de eucalipto submetidos a desbaste. [Evaluation of Weibull and hyperbolic functions for fitting data from eucalyptus stands under thinning]. Revista Árvore, Viçosa 34 (2): 305-311.

(32) Mabvurira D, Maltamo M, Kangas A (2002). Predicting and calibrating diameter distributions of Eucalyptus grandis (Hill) Maiden plantations in Zimbabwe. New Forests 23: 207-223.

(33) Mateus A, Tomé M (2011). Modelling the diameter distribution of eucalyptus plantations with Johnson’s SB probability density function: parameters recovery from a compatible system of equations to predict stand variables. Annals of Forest Science 68: 325-335.

(34) Maltamo M, Puumalainen J, Päivinen R (1995). Comparison of Beta and Weibull functions for modelling basal area diameter distributions in stands of Pinus sylvestris and Picea abies. Scandinavian Journal of Forest Research 10: 284-295.

(35) Merganic J, Sterba H (2006). Characterization of diameter distribution using the Weibull function: method of moments. European Journal of Forest Research 125: 427-439.

(36) Nabil MA (1989). Essai de synthèse sur la végétation et la phyto-écologie tunisiennes. I - Elément de botanique et de phyto-écologie. [Synthesis assay on the Tunisian vegetation and phytoecology. I - Botany and phytoecology element] (vol. 4-6). Faculté des Sciences, Tunis, Tunisia, pp. 247. Gscholar

(37) Nanang DM (1998). Suitability of the normal, log-normal and Weibull distributions for fitting diameter distributions of neem plantations in Northern Ghana. Forest Ecology and Management 103: 1-7.

(38) Nanos N, Montero G (2002). Spatial prediction of diameter distributions models. Forest Ecology and Management 161: 147-158.

(39) Nelson TC (1964). Diameter distribution and growth of loblolly pine. Forest Science 10: 105-115. Online | Gscholar

(40) Nogueira GS, Leite HG, Campos JCC, Carvalho AF, Souza AL (2005). Modelo de distribuição diamétrica para povoamentos de Eucalyptus sp. submetidos a desbaste [Diametric distribution model for thinned Eucalyptus sp. stands]. Revista Árvore, Viçosa 29 (4): 579-589.

(41) Palahi M, Pukkala T, Trasobares A (2006). Modelling the diameter distribution of Pinus sylvestris, Pinus nigra and Pinus halepensis forest stands in Catalonia using the truncated Weibull function. Forestry 79 (5): 553-562.

(42) Palahi M, Pukkala T, Blasco E, Trasobares A (2007). A comparison of beta, Johnson’s SB, Weibull and truncated Weibull functions for modelling diameter distribution of forest stands in Catalonia (north-east of Spain). European Journal of Forest Research 126 (4): 563-571.

(43) Reynolds MR, Burk TE, Huang WC (1988). Goodness-of-fit tests and model selection procedures for diameter distribution models. Forest Science 34: 373-399. Online | Gscholar

(44) Sevillano Marco E, Fernández-Manso A, Castedo-Dorado F (2009). Development and application of a growth model for Pinus radiata D. Don plantations El Bierzo (Spain). Investigation Agraria Sistemas Recursos Forestales 18 (1): 64-80. Gscholar

(45) Sghaier T, Palm R (2002). Répartition des arbres et des volumes par classes de grosseurs dans les peuplements de pin d’Alep (Pinus halepensis Mill.) en Tunisie. [Distribution of trees and volumes by diameter classes of the Aleppo pine (Pinus halepensis Mill.) stands in Tunisia]. Annals of Forest Science 59: 293-300.

(46) Sghaier T, Tomé M, Tomé J, Sánchez-González M, Cañellas I, Calama R (2013). Distance-independent individual tree diameter-increment model for Thuya (Tetraclinis articulata (VAHL.) MAST.) stands in Tunisia. Forest Systems 22 (3): 433-441.

(47) Sghaier T, Sánchez-González M, Garchi S, Ammari Y, Cañellas I, Calama R (2015). Developing a stand based growth and yield model for thuya (Tetraclinis articulata (Vahl) Mast) in Tunisia. iForest - Biogeosciences and Forestry 9 (1): 79-88.

(48) Tenbergen B, Günster A, Schreiber KF (1995). Harvesting runoff: the minicatchment technique - An alternative to irrigated tree plantations in semiarid regions. Ambio 24 (2): 72-76. Gscholar

(49) Zhang L, Packard KC, Liu C (2003). A comparison of estimation methods for fitting Weibull and Johnson’s SB distributions to mixed spruce-fir stands in northeastern North America. Canadian Journal of Forest Research 33: 1340-1347.