Modeling Diameter Distributions with Six Probability Density Functions in Pinus halepensis Mill. Plantations Using Low-Density Airborne Laser Scanning Data in Arag ó n (Northeast Spain)

: The diameter distributions of trees in 50 temporary sample plots (TSPs) established in Pinus halepensis Mill. stands were recovered from LiDAR metrics by using six probability density functions (PDFs): the Weibull (2P and 3P), Johnson’s SB, beta, generalized beta and gamma-2P functions. The parameters were recovered from the ﬁrst and the second moments of the distributions (mean and variance, respectively) by using parameter recovery models (PRM). Linear models were used to predict both moments from LiDAR data. In recovering the functions, the location parameters of the distributions were predetermined as the minimum diameter inventoried, and scale parameters were established as the maximum diameters predicted from LiDAR metrics. The Kolmogorov–Smirnov (KS) statistic ( D n ), number of acceptances by the KS test, the Cram é r von Misses ( W 2 ) statistic, bias and mean square error (MSE) were used to evaluate the goodness of ﬁts. The ﬁts for the six recovered functions were compared with the ﬁts to all measured data from 58 TSPs (LiDAR metrics could only be extracted from 50 of the plots). In the ﬁtting phase, the location parameters were ﬁxed at a suitable value determined according to the forestry literature (0.75 · d min ). The linear models used to recover the two moments of the distributions and the maximum diameters determined from LiDAR data were accurate, with R 2 values of 0.750, 0.724 and 0.873 for d g , d med and d max . Reasonable results were obtained with all six recovered functions. The goodness-of-ﬁt statistics indicated that the beta function was the most accurate, followed by the generalized beta function. The Weibull-3P function provided the poorest ﬁts and the Weibull-2P and Johnson’s SB also yielded poor ﬁts to the data.


Introduction
LiDAR (Laser Imaging Detection and Ranging) technology allows large areas to be scanned, thus generating continuous information about the entire space. In the last 20 years, LiDAR has been increasingly used in forest inventories at different scales [1,2], because of its ability to provide detailed three-dimensional information on the size and structure of forest cover [3]. The process provides data on the dimensions of the trees and on the structure of the forest cover, and other parameters such as canopy fuel attributes [4]. The variables most closely related to LiDAR data are canopy cover and tree height, which mainly depend on the vertical distribution of the layers forming the tree cover [5,6].
The main advantage of using field data that will then be processed using LiDAR techniques is that a much smaller sampling effort is required (with a good sampling design), regarding both the number of plots and the number of tree height measurements required. LiDAR technology, such as airborne laser scanning (ALS), provides a great from the data corresponding to 50 of these plots ( Figure 1). The plot size ranged from 225 to 625 m 2 , to achieve a minimum number of 30 trees per plot. Diameter at breast height (DBH at 1.3 m above the ground) of a total of 1685 trees was measured with callipers, to an accuracy of 0.1 cm. The following field and LiDAR variables were calculated from the inventory data: quadratic mean diameter, mean diameter, maximum diameter, number of trees per hectare, dominant height, basal area and total LiDAR return density within the plots (pulses·m −2 ). Summary statistics are shown in Table 1.
Remote Sens. 2021, 13, x FOR PEER REVIEW 3 of 17 other pines. Scant attention has therefore been given to Aleppo pine, even though at local scales it may be the only economic alternative available and may complement other nontimber harvests [22]. A total of 58 temporary sampling plots (TSPs) of P. halepensis Mill., measured in October 2017, were used for the present study, and LiDAR metrics were able to be extracted from the data corresponding to 50 of these plots ( Figure 1). The plot size ranged from 225 to 625 m 2 , to achieve a minimum number of 30 trees per plot. Diameter at breast height (DBH at 1.3 m above the ground) of a total of 1685 trees was measured with callipers, to an accuracy of 0.1 cm. The following field and LiDAR variables were calculated from the inventory data: quadratic mean diameter, mean diameter, maximum diameter, number of trees per hectare, dominant height, basal area and total LiDAR return density within the plots (pulses·m −2 ). Summary statistics are shown in Table 1.

Lidar Metrics
The LiDAR data covering 50 plots in Pinus halepensis stands were acquired in a national survey carried out for the PNOA (Plan Nacional de Ortofotografía de España) project between 16 October and 16 November 2016, under the direction of the Spanish Ministerio de Fomento (Dirección General del Instituto Geográfico Nacional and Centro Nacional de Información Geográfica). The data were obtained with a LEICA ALS80 sensor, at a pulse repetition rate of 45 kHz, scan frequency of 70 Hz and a maximum scan angle of ±50 • . A maximum of 4 returns per pulse were registered, reaching an average point density of 1.070 points·m −2 , with a theoretical laser pulse density required for the PNOA project of 0.5 first returns per square meter. Summary statistics of the LiDAR return density per square meter within the plots are shown in Table 1. The LiDAR data were also processed with FUSION software [23]. As reported by the provider (https://pnoa.ign.es/) accessed on 1 April 2020, the vertical accuracy of the LiDAR metrics, given by the RMSE, is ≤0.20 m. The set of metrics from the points laid above 2 m was extracted for each plot (Table 2). The three-parameter Weibull PDF has the following expression for a continuous random variable x [24]: where f(x) is the probability density of trees with diameter x, a represents the location, b the scale and c the shape. If a equals zero, the expression corresponds to the two-parameter Weibull function. The method of moments was used to fit all six distributions on the basis of the relationship between the parameters and the first and second moments of the diameter distribution (arithmetic mean diameter and variance). For the Weibull distribution, it is computed by the following expressions [10,18] Remote Sens. 2021, 13, 2307 5 of 17 where d is the arithmetic mean diameter of the distribution, σ 2 is the variance and Γ(·) is the gamma function. Equation (3) was resolved by a bisection iterative procedure [25].

The Beta Function
The general expression of the beta distribution for a random variable x is given by [26]: for the interval L ≤ x ≤ U, and x = 0 otherwise, where x is the diameter at breast height and is assumed to be continuous, f(x) is the density associated with diameter x, U and L are respectively the upper and lower limits of the beta distribution, c is the scaling factor, α and γ are respectively the first and the second exponents that determine the shape of the distribution and Γ(·) is the gamma function. The method of moments for the beta function [9,15,27] is computed by the following equations: where d is the arithmetic mean diameter of the distribution and s 2 is the variance.

The Generalized Beta Distribution (GBD)
The general expression of the generalized beta PDF for a random variable x is as follows [16]: for the interval (β 1 , β 1 + β 2 ), and 0 otherwise, where x is the diameter at breast height and is assumed to be continuous, f(x) is the density associated with diameter x, β 1 and β 2 are respectively the lower and upper limits of the distribution, β 3 and β 4 are exponents that determine the shape of the distribution, is the scaling factor of the function and Γ(·) is the gamma function. The method of moments for the GBD function is computed by the following expressions [16,28]: where x is the tree diameter, α 1 is the arithmetic mean diameter and α 2 is the variance.

The Johnson's SB Function
The model of the Johnson's SB PDF has the following expression for a continuous random variable x [29]: The model is characterized by the parameters ε (location), λ (scale), γ (asymmetry) and δ (kurtosis). The method of moments used in [11,30] is computed as follows: where d is the arithmetic mean of the plot diameters, Sd(x) is the modified standard (δ) deviation and σ x is the plot diameter standard deviation.

The Gamma Function
The model of the gamma PDF has the following expression for a continuous random variable x [14,31]: with x > γ, α > 0 and β > 0, where α is the shape parameter, β is an inverse scale parameter, γ is the location parameter (γ = 0 for the two-parameter gamma distribution) and Γ(·) is the gamma function. The method of moments for the gamma-2P function is computed by the following expressions [28]: where d is the mean diameter and s is the standard deviation. For fitting the distributions to the observed data (phase 1), a value of 0.75·d min was used for the location parameters, as [32] achieved good results with this value compared to other constraints. For the Weibull-2P and gamma-2P functions, the value of the location parameters was zero. Scale parameters and the upper limit of the four-parameter functions (Johnson's SB, beta and generalized beta) were established as the maximum diameter of the distributions (d max ) and the upper limit of the largest diameter class for the beta function [27,28] to improve parameters' convergence. The size of the diameter classes was assumed to equal 1 cm.

Goodness of Fit Evaluation
The consistency of the functions was evaluated using the Kolmogorov-Smirnov (KS) (D n ) and Cramér von Mises (W 2 ) statistics, bias, mean square error (MSE) and number of acceptances by the KS test. For a given cumulative distribution function, F(x), where sup x is the supremum of the set of distances, calculated as follows [33]: where the cumulative observed and estimated frequencies, i.e., F n (x i ) and F 0 (x j ), are compared. The recovered distributions were also assessed using the two-sample Kolmogorov-Smirnov test (KS) under the null hypothesis that the plot data originate from the recovered distribution. However, as the parameters were estimated empirically, the theoretical distribution for each plot is unknown, and the KS test should therefore be conducted using a Monte Carlo simulation [34]. The procedure is explained in further detail in [10,11].
The Cramér von Mises statistic (W 2 ) is a measure of the square of the distance between the empirical and the cumulative theoretical distribution [35]: whereF(x i ) is the cumulative theoretical distribution in diameter class i, and n is the number of diameter classes. The bias and the mean square error (MSE) were also used as goodness-of-fit measures and were expressed as follows: where Y i is the observed relative frequency of trees in each diameter class,Ŷ i is the theoretical value predicted by the model and N is the number of diameter classes. The bias and MSE were calculated for each fit as the mean relative frequency of trees.

Recovering the Parameters of the Distributions from LiDAR Metrics
The function parameters were recovered from the first two moments of the distributions (mean diameter and variance) by using parameter recovery models (PRMs). The variance σ 2 d is related to the quadratic mean diameter (d g ) and the mean diameter (d m ), as follows: Within the LiDAR data framework, the stand level attributes usually used for the recovery procedure (e.g., t, N, H, S) can be replaced by LiDAR metrics as explanatory variables to estimate d g and d m [10]. However, in the four-parameter functions (Johnson's SB, beta and generalized beta), the range and the upper limit of the distributions can also be estimated. In this case, both values were used as the maximum diameter (d max ), which was then related to LiDAR metrics. Location parameters of the recovered functions were predetermined as the minimum diameter inventory (7.5 cm) for the Weibull-3P, Johnson's SB, beta and generalized beta functions, and as zero for the Weibull-2P and gamma-2P functions.
We used linear models to establish the empirical relationships between d max , d g and LiDAR metrics: where d max and d g are the dependent variables, X 1 , X 2 , . . . , X m represents the independent, potentially explanatory variables related to the LiDAR-derived height distribution and canopy closure, α 0 , α 1 , . . . α n and β 0 , β 1 , . . . β n are the parameters to be estimated in the fitting process and ε is the additive error term, which is assumed to be independent and normally and identically distributed, with zero mean. For a given stand, d m is always smaller than or equal to d g , and we therefore used the following model expression to take this restriction into account [9][10][11]36]: where X 1 , X 2 , . . . , X m are the potential explanatory variables related to the LiDAR-derived height distribution and canopy closure, δ 0 ,δ 1 , . . . δ m are the parameters to be estimated in the fitting process and ε is the additive error term. Equation (30) was linearized by applying a natural logarithmic transformation to facilitate selection of the independent variables. Finally, Equation (29) (when d g was the dependent variable) and Equation (30) were fitted simultaneously with "seemingly unrelated regression" (SUR) to prevent cross-correlation between error components. Goodness of fits were evaluated with the coefficient of determination (R 2 ) and the root mean square error (RMSE).

Results
The mean, standard deviation, minimum and maximum values of the parameters estimated in step 1 (fitting the distributions to observed data) and step 2 (recovering the parameters of the distributions from LiDAR metrics) are shown in Table 3.
The mean values of the parameters were reasonable in both cases, with no large differences in the values produced by each method. Table 4 shows the mean value of the statistics used to test the goodness-of-fits to the observed distributions: Kolmogorov-Smirnov (D n ), Cramér von Mises (W 2 ), bias and mean squared error (MSE).
In the fitting phase, the lowest value of D n was yielded by the beta function (0.138560), followed by the Weibull-2P function (0.146922) and the Weibull-3P function (0.150761). The highest value was yielded by the gamma-2P function (0.164750), followed by the Johnson's SB function (0.162578) and the generalized beta function (0.153656). Regarding the Cramér von Mises statistic (W 2 ), the lowest value corresponded to the Johnson's SB function (0.043088) and the highest to the Weibull-2P function (0.081863). The order of the functions in terms of W 2 was Johnson's SB < Weibull-3P < generalized beta < gamma-2P < beta < Weibull-2P. The smallest MSE value corresponded to the beta function (0.001851), followed by the generalized beta function (0.001879). The highest MSE value was yielded by the Weibull-2P function (0.002009). The order for the six distributions studied was beta < generalized beta < Johnson's SB < Weibull-3P < gamma-2P < Weibull-2P. The bias may be less important for comparison of the results because errors with different signs can compensate each other in the total mean value. Thus, considering the values of D n , W 2 and MSE for the six distributions studied, the beta function provided the most accurate fits to the observed data and the Weibull-2P and the Gamma-2P functions yielded the poorest fits.
The parameter estimates and goodness-of-fit statistics (R 2 , RMSE and % RMSE) of the simultaneous fitting of Equations (29) and (30) used to estimate the mean diameter (d m ) and the quadratic mean diameter (d g ) from LiDAR data and for Equation (28) used to estimate d max are shown in Table 5. Table 3. Descriptive statistics for the parameters of the functions in the fitting and recovery steps.

Function
Step   The maximum diameter of the distributions (d max ) was considered the scale parameter for the Weibull-3P, Johnson's SB and generalized beta functions and for establishing the upper limit of the beta function. The maximum diameter was predicted by the following independent variables: mean absolute deviation for height (LH_AAD) and 95% height percentile (LH_P95). The linear model yielded an accuracy of R 2 = 0.873 and RMSE = 2.758. Models [29,30] used for simultaneous fitting of quadratic mean diameter (d g ) and the mean diameter (d med ) yielded R 2 values of 0.750 and 0.724 and RMSE values of 2.542 and 2.549, respectively. Good results were obtained for these key variables in terms of recovering the diameter distributions with the six functions.
For recovering the distributions from LiDAR data (Table 6), the smallest value of D n was obtained by the beta function (0.254078), as in the fitting step, and then by the generalized beta function (0.264840) and the gamma-2P function (0.272107). The highest value was yielded by the Weibull-3P function (0.286170), followed by the Weibull-2P function (0.282500) and the Johnson's SB function (0.273945). The lowest value of the Cramér von Mises statistic (W 2 ) was also yielded by the beta function (0.381916) and the highest by the Weibull-2P function (0.498277). The order of the distributions for the W 2 value was beta < generalized beta < gamma-2P < Johnson's SB < Weibull-3P < Weibull-2P. As for D n and W 2 , the smallest value for the MSE corresponded to the beta distribution (0.002649), followed by the generalized beta function (0.002711). The highest value of the MSE corresponded to the gamma-2P function (0.002792). The order in terms of MSE for the six distributions recovered was beta < generalized beta < Weibull-2P < Weibull-3P < gamma-2P < Johnson's SB. As for the other statistics, bias was also higher for recovery of the distributions from LiDAR metrics than for fitting to the observed data. Results for the KS test are also consistent with the other statistics, with the following order for plots accepted: beta function (35 plots: 70%), gamma-2P function (34 plots: 68%), generalized beta function (33 plots: 66%), Johnson's SB function (32 plots: 64%), Weibull-2P function (28 plots: 56%) and Weibull-3P function (26 plots: 52%). Thus, considering the values of D n , W 2 , MSE and KS test for the six distributions studied, the beta function was the most accurate for recovering the parameters from LiDAR data, followed by the generalized beta function. The Weibull-3P, Weibull-2P and the Johnson's SB functions yielded poorer results.
The mean values of the MSE in each diameter class for the fits to the observed data and for recovering from LiDAR data for the best (beta) and the poorest (Johnson's SB) functions are shown in Figure 2a, while those for the generalized beta vs. Weibull-3P and for the gamma-2P vs. Weibull-2P functions are also shown in Figure 2b,c. Figure 2a shows that the MSE for the fits to the observed data was almost the same for both functions (beta and Johnson's SB) over the diameter range studied, decreasing when the diameters increased. However, significant differences were observed for recovery of the functions. The recovered distributions yielded much larger values of MSE, up to 25 cm, after which both values (fitted and recovered) were more similar. The recovered Johnson's SB yielded the largest values up to 19 cm, while the beta function was less accurate at between 19 and 24 cm. The difference between fitting and recovery of the smallest diameter classes increased with the number of data points.
The lowest value of the Cramér von Mises statistic (W 2 ) was also yielded by the beta function (0.381916) and the highest by the Weibull-2P function (0.498277). The order of the distributions for the W 2 value was beta < generalized beta < gamma-2P < Johnson's SB < Weibull-3P < Weibull-2P. As for Dn and W 2 , the smallest value for the MSE corresponded to the beta distribution (0.002649), followed by the generalized beta function (0.002711). The highest value of the MSE corresponded to the gamma-2P function (0.002792). The order in terms of MSE for the six distributions recovered was beta < generalized beta < Weibull-2P < Weibull-3P < gamma-2P < Johnson's SB. As for the other statistics, bias was also higher for recovery of the distributions from LiDAR metrics than for fitting to the observed data. Results for the KS test are also consistent with the other statistics, with the following order for plots accepted: beta function (35 plots: 70%), gamma-2P function (34 plots: 68%), generalized beta function (33 plots: 66%), Johnson's SB function (32 plots: 64%), Weibull-2P function (28 plots: 56%) and Weibull-3P function (26 plots: 52%). Thus, considering the values of Dn, W 2 , MSE and KS test for the six distributions studied, the beta function was the most accurate for recovering the parameters from LiDAR data, followed by the generalized beta function. The Weibull-3P, Weibull-2P and the Johnson's SB functions yielded poorer results.
The mean values of the MSE in each diameter class for the fits to the observed data and for recovering from LiDAR data for the best (beta) and the poorest (Johnson's SB) functions are shown in Figure 2a, while those for the generalized beta vs. Weibull-3P and for the gamma-2P vs. Weibull-2P functions are also shown in Figure 2b,c. Figure 2a shows that the MSE for the fits to the observed data was almost the same for both functions (beta and Johnson's SB) over the diameter range studied, decreasing when the diameters increased. However, significant differences were observed for recovery of the functions. The recovered distributions yielded much larger values of MSE, up to 25 cm, after which both values (fitted and recovered) were more similar. The recovered Johnson's SB yielded the largest values up to 19 cm, while the beta function was less accurate at between 19 and 24 cm. The difference between fitting and recovery of the smallest diameter classes increased with the number of data points.   Graphical assessment of the observed, described (fitted) and recovered distributions for six plots fits are shown in Figure 3. Pairwise comparisons of the functions were made: beta vs. Johnson's SB (Figure 3a), generalized beta vs. Weibull-3P (Figure 3b) and gamma-2P vs. Weibull-2P (Figure 3c). The predictions were reasonable and followed the observed diameter distributions, especially for plots 1, 3, 4 and 6.  Figure 2b shows the same trends for the range studied, with similar behavior obtained for fitting and recovery of both the generalized and Weibull-3P distributions. Figure 2c shows similar results produced by the gamma-2P and Weibull-3P functions.
Graphical assessment of the observed, described (fitted) and recovered distributions for six plots fits are shown in Figure 3. Pairwise comparisons of the functions were made: beta vs. Johnson's SB (Figure 3a), generalized beta vs. Weibull-3P (Figure 3b) and gamma-2P vs. Weibull-2P (Figure 3c). The predictions were reasonable and followed the observed diameter distributions, especially for plots 1, 3, 4 and 6. emote Sens. 2021, 13, x FOR PEER REVIEW 13 of 17

Discussion
We used a total of 58 TSPs from P. halepensis Mill. plantations for the present study, and LiDAR metrics were able be extracted from 50 of these plots. The number of plots is greater than in other studies of diameter distributions recovered from LiDAR data in the Northwest Iberian Peninsula [10,11]. These studies recovered the Weibull-2P and the Johnson's SB distributions from LiDAR metrics. The Weibull function (2P and 3P) is the most commonly used distribution in this type of study [20,37,38]. However, in the present study, we used, for the first time within the framework of LiDAR-based research, three functions commonly used to describe and predict diameter distributions in forest stands, i.e., the beta [26], generalized beta [16] and gamma-2P [14] functions. In some cases, the results obtained were better than those obtained with the more usual Weibull (2P and 3P) and Johnson's SB functions, for both fitting to the observed data and recovering with LiDAR metrics.
For the fitting phase, the results obtained for the main statistics used (D n , W 2 and MSE) are similar to those obtained in [28] for Eucalyptus globulus stands in NW Spain. However, more accurate results were obtained in the same study for Pinus radiata stands by using the same functions. In a previous study, similar results were also obtained for Pinus sylvestris and better than those reported for Pinus pinaster stands in NW Spain [32].
The models used to recover the parameters of the distributions from moments of the distributions were accurate. Explanatory variables were LH_P90 (height percentile of 90%) for d g , and LH_MAD_MEDIAN (median of the absolute deviations from the overall height median) for d med . In other studies, height percentiles were also used as independent variables for estimating d g (75% percentile and number of LiDAR last returns above a height of 1 m) and for d med (1% percentile) in the models obtained in [10] for recovering the Weibull-2P distribution in 25 plantation plots of Pinus radiata in Northwest Spain, with R 2 for the d g and d med models of 0.80 and 0.77, respectively. For Pinus halepensis, we obtained similar R 2 values (0.75 and 0.72 for d g and d med , respectively). The RMSE values (2.54 cm for d g and 2.55 cm for d med ) are consistent with the values reported in international forestry literature. For example, for German forests dominated by Picea abies (L.) Karst., the authors of [39] used data from a 0.44 pulse m −2 LiDAR flight and reported an RMSE of 2.44 cm for d med , while the authors of [40] studied a broad range of forest types (coniferous and hardwoods) and conditions across Ontario by using an artificially reduced LiDAR database of 0.5 pulses m −2 , reporting RMSE values ranging from 0.76 to 4.3 cm for d g . The authors of [10] reported RMSE values of 3.42 for d g and 3.63 for d med , while the authors of [11] used exponential models instead of linear models to estimate the moments of the Johnson's SB and the Weibull-2P functions, obtaining an R 2 of 0.82 and 0.86 for the d med of Pinus radiata and Eucalyptus globulus respectively, and 0.84 and 0.89 for d g of the same species. We compared the use of linear and exponential models for obtaining both variables and found that the results were similar, with the exponential model even including more independent variables. Thus, the linear models were considered more suitable.
Use of the four-parameter Johnson's SB, beta and generalized beta distributions requires knowledge about the location parameter of the three distributions, the scale for the Johnson's SB and the upper limit of the distributions for the beta and generalized beta functions. The three-parameter Weibull-3P requires knowledge about the location parameter, while in the Weibull-2P and gamma-2P functions, the value of the location parameter is zero. In functions with a location parameter, the value considered in the recovery phase with LiDAR metrics was the minimum diameter inventory (7.5 cm). This avoided having to model this parameter with LiDAR metrics, which can produce inaccurate results if stand variables are used in the model [9]. Thus, in all cases, the location parameter was considered equal to 7.5 cm for recovering the distributions from LiDAR data, while for fitting to the observed data, it was assumed to be 0.75·d min , producing a similar mean value considering all plots of 7.18 cm ( Table 3). The authors of [11] considered the location parameter for the Johnson's SB to equal zero, thus avoiding modeling it with LiDAR metrics, while the scale parameter considered was also the maximum diameter (d max ).
The model obtained for d max included, as LiDAR explanatory variables, LH_AAD (mean absolute deviation for height) and LH_P95 (95% height percentile), with R 2 = 0.87 and % RMSE = 10.56. The authors of [11] obtained similar values for Pinus radiata (R 2 = 0.93 and % RMSE = 8%) and Eucalyptus globulus (R 2 = 0.83 and % RMSE = 12%). The results for this variable modeled from LiDAR data are more accurate than those obtained in [27] and [9] with stand variables. The statistics obtained for recovery of the six functions were also reasonable. For example, the authors of [32] reported D n values of 0.1924 for the Weibull-2P fitted by Maximum Likelihood, 0.2285 for the Johnson's SB fitted by conditional maximum likelihood (CML) and 0.1812 for the beta fitted by moments to observed distributions of Pinus pinaster in Northwest Spain. The authors of [33] obtained a D n value of 0.193 in the fits by maximum likelihood of the Weibull-3P to observed data of Pinus taeda plantations in the USA. The value obtained for the recovered beta from LiDAR metrics (0.2541) is close to these values. The number of KS acceptances are consistent with those found in previous studies. The percentage of acceptance (between 52% and 70%) was similar for Eucalyptus globulus in northwest Portugal [11] and higher than in Pinus radiata stands in northwest Spain [10,11].

Conclusions
Scant attention has been given to Aleppo pine plantations in Spain. However, we recovered the parameters of six probability density functions, i.e., the Weibull (2P and 3P), Johnson's SB, beta, generalized beta and gamma-2P functions, from LiDAR metrics in 50 Pinus halepensis stands in Aragón (northeast Spain). The beta, generalized beta and gamma-2P functions are of novel application in the field of LiDAR-based research. The results obtained were reasonable compared with previous studies and they showed that the beta and generalized beta functions were more accurate than the Johnson's SB and Weibull (2P and 3P) functions, which are more commonly used in this type of study. None of the distributions predicted anomalous results.