Symmetric and Asymmetric Diffusions through Age-Varying Mixed-Species Stand Parameters

Abstract

(1) Background: This paper deals with unevenly aged, whole-stand models from mixed-effect parameters diffusion processes and Voronoi diagram points of view and concentrates on the mixed-species stands in Lithuania. We focus on the Voronoi diagram of potentially available areas to tree positions as the measure of the competition effect of individual trees and the tree diameter at breast height to relate their evolution through time. (2) Methods: We consider a bivariate hybrid mixed-effect parameters stochastic differential equation for the parameterization of the diameter and available polygon area at age to ensure a proper description of the link between them during the age (time) span of a forest stand. In this study, the Voronoi diagram was used as a mathematical tool for the quantitative characterization of inter-tree competition. (3) Results: The newly derived model considers bivariate correlated observations, tree diameter, and polygon area arising from a particular stand and enables defining equations for calculating diameter, polygon-area, and stand-density predictions and forecasts. (4) Conclusions: From a statistical point of view, the newly developed models produced acceptable statistical measures of predictions and forecasts. All the results were implemented in the Maple computer algebra system.

1. Introduction

Forest statisticians often need to address complex issues in important whole-stand, unevenly aged growth models with high levels of uncertainty that affect individual trees. Providing scientific evidence for effective decision processes in forest areas is a key issue, and stochastic calculus is an essential tool [1].

In the past few decades, mathematical models based on diffusion processes were fruitfully applied to describe stochasticity phenomena belonging to even extremely different disciplinary fields that range in scale from human population in biology [2] to cryptocurrency in finance [3] and from infectious disease in medicine [4] to networks in neuroscience [5]. Mixed-effect parameters diffusion processes provide a convenient tool to account for differences between several experiments or several subjects [6].

One of the important aims of modern forest regeneration is the investigation of the mapped tree community distribution in a stand area and its evolution through time. The distribution of individual trees depends both on the initial conditions present in a forest stand and the successive development of trees through changes in size components and mortality events. At present, results are large datasets of observations and statistical analysis acquired by using competition indices to both direct and indirect stand measurements [7]. To operate with these data, it is necessary to investigate information in a convenient form for the mathematical modeling of individual-based tree growth and yield. Independently of the type of an individual tree or whole-stand stochastic modelling (having a random variable), the consideration of a diffusion process has particular challenges. Diffusion processes allow for us to compute both the solution of the corresponding stochastic differential equation in a form of probability density function and their statistical moments, mainly mean, quantiles, and variance. The mechanism of the evolution of the competition among living trees in a forest stand can be successfully described by using a relatively young mathematical tool—continuous stochastic differential equations (SDEs). SDEs serve the purpose of analysis in a wide variety of generic growth models with different resolutions. The most advanced statistical models of this kind were formulated in terms of multivariate diffusion processes [8].

The process by which tree seedlings start to grow successfully in a new place is spatially multivariate, with numerous unknown variables. A potentially available area for an individual seedling embodies a variable representing the competitive strength of individual tree–tree interactions and varies through time. Several case studies were conducted to study the relationships between the potentially available area of a tree and its main biometrical attributes and competition indices [9]. The competition index is widely applied in plant ecology, but it has a smaller number of applications in individual-tree or whole-stand modeling. From a mathematical point of view, the Voronoi diagram [10] can be regarded as the best solution for neighboring effects on tree-size growth. To illustrate this study, we use the Voronoi diagram method, which is used by the individual-based tree-growth model. The Voronoi diagram formalizes a detailed structure of the position, size, and shape of a potentially available area of individual trees with respect to the number of trees per unit area and nearness of their contiguous neighbors. Moreover, Voronoi polygons reflect the local variation of the number of trees in a given location. A variable outlined by the value of the area of the Voronoi polygon is used as a descriptive parameter of spatial arrangement and as a predictor of stand density models.

Almost all the published research on natural mechanisms of maintaining tree-species diversity focuses on the principles of spatial measures and competition-index construction and their biological and mathematical reasoning [9]. Previous relationships on the key competition index include only local neighbors and are static. Unfortunately, it is insufficient, as each competition index systematically changes through time, which is of great importance and highly correlated with tree-size variables. Therefore, this paper focuses on the Voronoi polygon area of a tree and its dynamics. The dynamics of the polygon area is described by a diffusion process that accounts for the competition effect of contiguous neighbors and its correlation with tree-size or -stand variables.

This paper deals with whole-stand models from a diffusion process of view and concentrates on mixed-species stands in Lithuania. The innovation of this paper lies in the following: (1) the Voronoi diagram of potentially available areas to tree positions as the measure of the competition effect of individual trees, and its evolution through time is described using a Gompertz-type diffusion process; (2) the link of stand density with potentially available area and tree diameter is analyzed; (3) the mean and quantile attribute equations of the potentially available area, tree diameter, and stand density are described; and (4) the significance of differences in the distributions of different tree species is assessed. This paper carries out further research from the following directions: the first aspect analyzes the stochastic fixed-effect parameters process and then the stochastic mixed-effect parameters process; the third aspect illustrates the computational properties of the derived models.

The study was conducted in naturally and artificially regenerated areas by using data collected from over 58,000 trees positioned in a network of 50 permanent, rectangular sampling plots. As significant variables in the present study, we used such assessed attributes as species, the location of the individuals (Cartesian coordinates x, y), tree diameter at 1.3 m above ground level, and age. In order to determine the individual area of the i-th (i = 1, …, n) tree, Voronoi polygon area P_i was calculated.

By utilizing the additive and multiplicative noise systems, we define the hybrid bivariate Vasicek–Gompertz-type SDE. The inclusion of random effects in the SDE produced a more accurate model to account for the error between model predictions and observed dataset. Model building in forestry was traditionally formalized by empirical models supporting the decisions of forest managers [11]. In addition, large published growth and yield systems were described by regression analysis, developed on a statistical basis that actually sought only to highlight the importance of a particular response variable in predicting the future, remaining static and deterministic. The stochastic differential equation approach opened the way to a stochastic dynamical formulation of the individual-tree or whole-stand growth process that appeared much later in statistical forestry [12,13,14].

2. Materials and Methods

2.1. Voronoi Diagram

Considering data in the form of a given set of n points (x_i; y_i), i = 1, …, n within a planar region A, we can assign an area to (x_i; y_i) consisting of that part of A that is closer to (x_i; y_i) than to any other point (x_j; y_j). A Voronoi diagram involves the partitioning of a planar region A into regions or polygons A_i on the basis of the distance to a specified discrete set of points [15]. The Voronoi diagram is a mathematical tool for the quantitative characterization of natural phenomena, such as inter-tree competition. In these polygons, trees with which the area neighbors are in direct competition for available light and nutrients are shown [16]. Each area A_i is a convex polygonal region. Voronoi polygons A_i vary with increasing time, and the number of trees decreases. Considering the dynamic of the Voronoi diagram allowed for both the quantification of its temporal dynamics and the characterization of the number of trees in a stand and its dynamics through time.

A generic definition of Voronoi diagrams was mathematically first introduced by Aurenhammer [17]. Voronoi cell A_i, associated with point (x_j; y_i), is the set of all points in A, of which distance to point (x_j; y_i) is not greater than their distance to other point (x_j; y_j), where j is any index different from i. If we assume that d((x; y),B) = inf{d((x; y),(a₁; a₂))|(a₁; a₂) ∈B} denotes the distance between point (x; y) and subset B, then:

$$A_i=\left\{(x;y)\in A|d\left((x;y), A_i\right)\le d\left((x;y), A_j\right) for all j \ne i\right\}.$$

(1)

The Voronoi diagram is simply the corpus of polygons A_i, i = 1, …, n. In the past decade, mathematical algorithms and the rapid development of computational capacities provided a new chance to develop Voronoi diagram-based applications. In this study, we discuss mathematical models and their analyses based on areas of generated Voronoi polygons. To illustrate our work, we computed the Voronoi polygon of each individual tree in a plot, which was remeasured four times and is visualized in Figure 1, which shows that the areas of the Voronoi cells change through age (time).

2.2. Bivariate SDEs of Diameter and Polygon Area

In forestry studies, most tree attributes are considered to be functions of the tree diameter at breast height D [18]. For explaining the dynamics of tree-size variables, the incorporation of the distance of neighboring trees into a model improves predictions [19]. Moreover, appropriate evaluation of individual-tree size-variable is a fundamental requirement for the analysis of stand variables, such as the number of trees and volume per hectare. This study considers a stochastic approach for the parameterization of diameter D and polygon area P at age to ensure the proper descriptions of their link during the age (time) span of a forest stand. In this section, we state the stochastic mixed-effect bivariate (diameter D and polygon area P) model and derive the corresponding bivariate probability density function. We also introduce the procedure for random-effect calibration. The Vasicek–Gompertz-type bivariate hybrid stochastic mixed-effect parameters model was used to parameterize diameter and polygon area at age-discrete data. The univariate Gompertz-type model was applied to the analysis of stem volume and tree-stem taper from different tree species [20]. Below, stochastic vector $X^i\left(t\right)={\left(X_1^i\left(t\right),X_2^i\left(t\right)\right)}^T={\left(D^i\left(t\right),P^i\left(t\right)\right)}^T$, i = 1,…,M, $t\in \left[t_0;T\right]$, M is the number of individuals (plots), and T is a finite horizon, T < ∞. Hybrid bivariate Vasicek–Gompertz-type SDE is defined as:

$$d{X}^i\left(t\right)=A^i\left(X^i\left(t\right)\right)dt+D\left(X^i\left(t\right)\right)B^{\frac{1}{2}}\cdot d{W}^i\left(t\right)P\left(X^i\left(t_0\right)=x_0\right)=1, i=1, \dots , M,$$

(2)

$$A^i\left(x\right)={\left({\beta }_d\left({\alpha }_d+{\phi }_d^i-x_1\right),\left(\left({\alpha }_p+{\phi }_p^i\right)-{\beta }_{p}ln(x_2)\right)x_2\right)}^{T}B=\left(\begin{array}{cc}{\sigma }_{dd}& {\sigma }_{dp}\\ {\sigma }_{dp}& {\sigma }_{pp}\end{array}\right),$$

(3)

$$\left(D\left(x\right)B^{\frac{1}{2}}\right){\left(D\left(x\right)B^{\frac{1}{2}}\right)}^T=\left(\begin{array}{cc}1& 0\\ 0& x_2\end{array}\right)\left(\begin{array}{cc}{\sigma }_{dd}& {\sigma }_{dp}\\ {\sigma }_{dp}& {\sigma }_{pp}\end{array}\right)\left(\begin{array}{cc}1& 0\\ 0& x_2\end{array}\right)$$

This study focusses on initial distribution, defined by deterministic initial value $\left(X_1^i\left(t_0\right),X_2^i\left(t_0\right)\right)=(x_{10},x_{20})=\left(x_{10},\mathsf{\delta }+{\phi }_0^i\right)$, and δ is an unknown fixed-effect parameter to be estimated. SDE of form (2) consists of two parts. The deterministic part in the model is drift function $A^i\left(x\right)$. Random term $D\left(X^i\left(t\right)\right)B^{\frac{1}{2}}\cdot d{W}^i\left(t\right)$, which corresponds to the uncertain part of the model, is referred to as the system noise. In Equation (2), $W^i\left(t\right)={\left(W_1^i\left(t\right),W_2^i\left(t\right)\right)}^T$ represents the bivariate Brownian motion, of which the time derivative is white noise. Moreover, Brownian motion increments $d{W}^i\left(t\right)$, i = 1,…,M, are considered to be independent across all plots. Random effects ${\phi }_d^i$${\phi }_p^i$, and ${\phi }_0^i$ are independent and normally distributed random variables with zero mean and constant variances, respectively, ${\phi }_d^i~N\left(0;{\sigma }_d^2\right)$, ${\phi }_p^i~N\left(0;{\sigma }_p^2\right)$, and ${\phi }_0^i~N\left(0;{\sigma }_0^2\right)$. Fixed-effect parameters vector θ to be estimated is defined as:

$$\theta =\left({\alpha }_d, {\alpha }_p,{\beta }_d,{\beta }_p,{\sigma }_{dd}, {\sigma }_{pp},{\rho }_{dp},\delta ,{\sigma }_d,{\sigma }_p, {\sigma }_0\right){\sigma }_{dp}=\sqrt{{\sigma }_{dd}{\sigma }_{pp}}{\rho }_{dp}.$$

(4)

The solution of the Vasicek–Gompertz-type SDE (2) has bivariate normal-lognormal distribution $N_{1}L{N}_1\left({\mu }^i(t|t_0,x_0);\hspace{0.33em}\Sigma (t|t_0)\right)$ with mean vector ${\mu }^i(t|t_0,x_0)$ and variance–covariance matrix $\Sigma (t|t_0)$, defined as:

$$\begin{gathered} {\mu }^i\left(t|t_0,x_0\right)=\left(\begin{array}{c}{\mu }_d^i\left(t|t_0,x_{10}\right)\\ {\mu }_p^i(t|t_0,x_{20})\end{array}\right) \\ =\left(\begin{array}{c}({\alpha }_d+{\phi }_d^i)+\left(x_{10}-({\alpha }_d+{\phi }_d^i)\right)e^{-{\beta }_d(t-t_0)}\\ e^{-{\beta }_p(t-t_0)}ln(\delta +{\phi }_0^i)+\frac{1-e^{-{\beta }_p(t-t_0)}}{{\beta }_p}\left({\alpha }_p+{\phi }_p^i-\frac{{\sigma }_{pp}}{2}\right)\end{array}\right) \end{gathered}$$

(5)

$$\begin{gathered} \Sigma \left(t|t_0\right)=\left(\begin{array}{cc}{v}_{dd}\left(t|t_0\right)& v_{dp}(t|t_0)\\ v_{dp}(t|t_0)& v_{pp}(t|t_0)\end{array}\right) \\ =\left(\begin{array}{cc}\frac{1-e^{-2{\beta }_d\left(t-t_0\right)}}{2{\beta }_d}{\sigma }_{dd}& \frac{1-e^{-\left({\beta }_d+{\beta }_p\right)\left(t-t_0\right)}}{{\beta }_d+{\beta }_p}\sqrt{{\sigma }_{dd}{\sigma }_{pp}}{\rho }_{dp}\\ \frac{1-e^{-\left({\beta }_d+{\beta }_p\right)\left(t-t_0\right)}}{{\beta }_d+{\beta }_p}\sqrt{{\sigma }_{dd}{\sigma }_{pp}}{\rho }_{dp}& \frac{1-e^{-2{\beta }_p\left(t-t_0\right)}}{2{\beta }_p}{\sigma }_{pp}\end{array}\right). \end{gathered}$$

(6)

We can separately calculate the probability distribution of each random variable if we wish to restrict our attention to the value of just one, for example, diameter D or polygon area P. The marginal distribution of $X_1^i\left(t\right)|X_1^i\left(t_0\right)=x_{10}$ is normal $N_1\left({\mu }_d^i(t|t_0,x_{10}) ;\hspace{0.33em}{v}_{dd}\left(t|t_0\right)\right)$, and the marginal distribution of $X_2^i\left(t\right)|X_2^i\left(t_0\right)=x_{20}$ is lognormal $L{N}_1\left({\mu }_p^i(t|t_0,x_{20}) ;\hspace{0.33em}{v}_{pp}(t|t_0)\right)$ with means ${\mu }_d^i(t|t_0,x_{10})$ and ${\mu }_p^i(t|t_0,x_{20})$ and variances $v_{dd}\left(t|t_0\right)$ and $v_{pp}\left(t|t_0\right)$, respectively.

For the diameter dynamic, the mean, median, mode, qth quantile (0 < q < 1), and variance trends can be listed as:

$$m_d^i\left(t|t_0,x_{10}\right)=m{e}_d^i\left(t|t_0,x_{10}\right)=m{o}_d^i\left(t|t_0,x_{10}\right)={\mu }_d^i(t|t_0,x_{10}),$$

(7)

$$m{q}_d^i\left(t,q|t_0,x_{10}\right)={\Phi }_q^{-1}\left({\mu }_d^i\left(t|t_0,x_{10}\right);v_{dd}\left(t|t_0\right)\right),$$

(8)

$$w_d^i\left(t|t_0,x_{10}\right)=v_{dd}\left(t|t_0\right),$$

(9)

where ${\Phi }_q^{-1}\left(\cdot ;\cdot \right)$ is the inverse of the standard normal distribution function.

For the polygon area dynamic, the mean, median, mode, qth quantile (0 < q < 1), and variance trends can be listed as:

$$m_p^i\left(t|t_0,x_{20}\right)=exp\left({\mu }_p^i(t|t_0,x_{20})+\frac{1}{2}{v}_{pp}\left(t|t_0\right)\right),$$

(10)

$$m{e}_p^i\left(t|t_0,x_{20}\right)=exp\left({\mu }_p^i(t|t_0,x_{20})\right),$$

(11)

$$m{o}_p^i\left(t|t_0,x_{20}\right)=exp\left({\mu }_p^i\left(t|t_0,x_{20}\right)-v_{pp}\left(t|t_0\right)\right),$$

(12)

$$m{q}_p^i\left(t,q|t_0,x_{20}\right)=L{\Phi }_q^{-1}({\mu }_p^i(t|t_0,x_{20});v_{pp}\left(t|t_0\right)),$$

(13)

$$w_p^i\left(t|t_0,x_{10}\right)=exp\left(2{\mu }_p^i(t|t_0,x_{20})+v_{pp}\left(t|t_0\right)\right)\cdot \left(\exp \left(v_{pp}\left(t|t_0\right)\right)-1\right),$$

(14)

where $L{\Phi }_q^{-1}\left(\cdot ;\cdot \right)$ is the inverse of the standard normal distribution function.

Conditional distribution of $X_1^i\left(t\right)|X_1^i\left(t_0\right)=x_{10}$ at a given $\left(X_2^i(t)=x_2\right)$ is univariate normal $N_1\left({\eta }_d^i(t,x_2|t_0,x_0);\hspace{0.33em}{\lambda }_d(t|t_0)\right)$, and conditional distribution of $X_2^i\left(t\right)|X_2^i\left(t_0\right)=x_{20}$ at a given $\left(X_1^i\left(t\right)=x_1\right)$ is univariate lognormal $L{N}_1\left({\eta }_p^i(t,x_1|t_0,x_0);{\lambda }_p(t|t_0)\right)$, with means and variances given as:

$${\eta }_d^i\left(t,x_2|t_0,x_0\right)={\mu }_d^i(t|t_0,x_{10})+\frac{v_{dp}\left(t|t_0\right)}{v_{pp}\left(t|t_0\right)}\left(\ln \left(x_2\right)-{\mu }_p^i(t|t_0,x_{20}) \right)$$

(15)

$${\lambda }_d\left(t|t_0\right)=v_{dd}\left(t|t_0\right)-\frac{{\left(v_{dp}\left(t|t_0\right)\right)}^2}{v_{pp}\left(t|t_0\right)},$$

(16)

$${\eta }_p^i\left(t,x_1|t_0,x_0\right)={\mu }_p^i(t|t_0,x_{20})+\frac{v_{dp}\left(t|t_0\right)}{v_{dd}\left(t|t_0\right)}\left(x_1-{\mu }_d^i(t|t_0,x_{10})\right),$$

(17)

$${\lambda }_p\left(t|t_0\right)=v_{pp}\left(t|t_0\right)-\frac{{\left(v_{dp}\left(t|t_0\right)\right)}^2}{v_{dd}\left(t|t_0\right)}.$$

(18)

For the diameter growth model, the conditional mean, median, mode, qth quantile (0 < q < 1), and variance trends can be listed as:

$$m{c}_d^i\left(t,x_2|t_0,x_0\right)=me{c}_d^i\left(t,x_2|t_0,x_0\right)=mo{c}_d^i\left(t,x_2|t_0,x_0\right)={\eta }_d^i\left(t,x_2|t_0,x_0\right),$$

(19)

$$mq{c}_d^i\left(t,x_2,q|t_0,x_0\right)={\Phi }_q^{-1}({\eta }_d^i(t,x_2|t_0,x_0);\hspace{0.33em}{\lambda }_d\left(t|t_0\right)),$$

(20)

$$w{c}_d^i\left(t|t_0\right)={\lambda }_d\left(t|t_0\right).$$

(21)

For the polygon area growth model, the conditional mean, median, mode, qth quantile (0 < q < 1), and variance trends can be listed as:

$$m{c}_p^i\left(t,x_1|t_0,x_0\right)=exp\left({\eta }_p^i\left(t,x_1|t_0,x_0\right)+\frac{1}{2}{\lambda }_p\left(t|t_0\right)\right),$$

(22)

$$me{c}_p^i\left(t,x_1|t_0,x_0\right)=exp\left({\eta }_p^i\left(t,x_1|t_0,x_0\right)\right),$$

(23)

$$mo{c}_p^i\left(t,x_1|t_0,x_0\right)=exp\left({\eta }_p^i\left(t,x_1|t_0,x_0\right)-{\lambda }_p\left(t|t_0\right)\right),$$

(24)

$$mq{c}_p^i\left(t,x_1,q|t_0,x_0\right)={\Phi }_q^{-1}({\eta }_p^i(t,x_1|t_0,x_0);\hspace{0.33em}{\lambda }_p\left(t|t_0\right)),$$

(25)

$$w{c}_p^i\left(t,x_1|t_0,x_0\right)=exp\left(2{\eta }_p^i\left(t,x_1|t_0,x_0\right)+{\lambda }_p\left(t|t_0\right)\right)\cdot \left(exp\left({\lambda }_p\left(t|t_0\right)\right)-1\right).$$

(26)

2.3. Data

The field-study area is located in the municipality of Kazlų Rūda in Lithuania. A major part of the Kazlų Rūda municipality is located in the fertile Užnemunė lowland and is among the most wooded areas in Lithuania, with about 59.4% of the territory covered by large forests. The specific allocation comprises the area covered by stands of pine (Pinus sylvestris), 63.8%; spruce (Picea abies), 30.2%; silver birch (Betula pendula Roth and Betula pubescens Ehrh.), 5.8%; and others, 0.2%. All data were collected during fieldwork in 1983–2019 across the Kazlų Rūda forests (latitude, 54°44′54.222″ N; longitude, 23°29′27.7944″ E; altitude, 68 m). Mean temperatures vary from −16.4 °C in winter to 22 °C in summer. Precipitation is distributed throughout the year, although predominantly in summer; the average is approximately 680 mm a year. During the 1983–1987 period, 50 permanent experimental plots were established in the Kazlų Rūda forests in Lithuania. According to regeneration mode, the 50 field-sample plots vary between those naturally and artificially regenerated and spread in pure or mixed-species stands. Each sample plot consisted of about 0.16–0.72 ha and was remeasured several times, from 1 until 6 (see Figure 1, showing 4 cycles) at 2- to 36-year intervals. The attributes recorded for each tree in the considered plots were the tree species, age, location of the sample trees (planar coordinate position x and y), and diameter at breast height. The age of the i-th tree (i varied from all the trees until the 10th) in the first measurement was recorded by counting its growth rings on the increment core (for even-aged stands, from entries in documents), and the ages of the remaining trees were obtained from the arithmetic mean. The accuracy of planar coordinate position was 1 dcm, and diameter measurements were performed with approximately 1 mm accuracy. The 50 field-sample plot dataset was randomly divided into estimation and validation datasets. The estimation dataset consisted of measurements from 40 plots, and the validation dataset compounded the remaining measurements from 10 plots.

Table 1. Tree characteristics of sample plots (0.16–0.72 ha) selected for the study.

Species	Data	Number of Trees	Min	Max	Mean	St. Dev.	Number of Trees	Min	Max	Mean	St. Dev.
Estimation							Validation
Pine	t (year)	28,982	5.0	172.0	49.14	22.16	6997	7.0	197.0	67.90	19.06
	d (cm)	28,982	0.5	61.0	17.87	9.17	6997	3.0	59.2	24.35	9.02
	p (m2)	28,982	0.09	120.21	9.24	7.56	6997	0.44	84.84	13.02	8.54
Spruce	t (year)	11,493	12.0	207.0	59.63	21.57	7001	7.0	191.0	67.94	21.01
	d (cm)	11,493	0.20	62.0	11.68	7.34	7001	3.0	61.8	13.24	8.51
	p (m2)	11,493	0.11	160.24	9.17	8.32	7001	0.26	77.76	9.46	7.65
Birch	t (year)	2880	5.0	107.32	47.75	17.88	663	16.0	129.73	60.29	15.87
	d (cm)	2880	0.90	45.40	14.46	8.13	663	3.0	50.0	19.91	9.30
	p (m2)	2880	0.33	173.82	8.70	7.01	663	0.89	51.93	9.88	8.05
	t (year)	43,410	5.0	207.0	51.84	22.27	14,711	7.0	197.0	67.55	19.94
All	d (cm)	43,410	0.20	62.0	16.0	9.07	14,711	3.0	61.80	18.87	10.35
	p (m2)	43,410	0.09	173.82	9.18	7.73	14,711	0.26	84.84	11.16	8.29

shows the summary statistics of the field-sample data.

3. Results

3.1. Parameter-Estimating Results

The key research problem in this study is the latest developments of accurate and computationally optimal parameter-estimation procedures based on the approximated maximum-likelihood technique, which cannot be implemented in the presence of a closed-form expression for the mixed-effect parameters SDE [21]. Our developed maximum-likelihood estimation technique relies on the fact that the conditional bivariate probability density function has an exact form. Therefore, the likelihood function maximization technique with respect to parameter vector θ for a given set of discretely observed diameter and polygon area data is defined by a two-step procedure.

To evaluate the proposed Vasicek–Gompertz-type mixed-effect parameters SDE (2), an approximated log-likelihood technique for estimating the parameters was set up on the basis of the discrete observations from estimation dataset $\left\{\left(d_1^i,p_1^i\right),\hspace{0.33em}\left(d_2^i,p_2^i\right),\hspace{0.33em}\dots ,\hspace{0.33em}\left(d_{n_i}^i,p_{n_i}^i\right)\right\}$ on a fixed time interval $\left\{t_1^i,t_2^i,\dots ,t_{n_i}^i\right\},$ i = 1,…,M. The randomly selected 40 samples were used to fit the SDE model defined by Equation (2), and the results of the estimating parameters are summarized in

Table 2. Vasicek–Gompertz-type system for diameter and polygon area: parameter estimates.

Species	${\mathit{\alpha }}_{\mathit{d}}$	${\mathit{\beta }}_{\mathit{d}}$	${\mathit{\alpha }}_{\mathit{p}}$	${\mathit{\beta }}_{\mathit{p}}$	${\mathit{\sigma }}_{\mathit{d}\mathit{d}}$	${\mathit{\sigma }}_{\mathit{p}\mathit{p}}$	${\mathit{\rho }}_{\mathit{d}\mathit{p}}$	δ	σd	σp	σ0
All	48.3358	0.0086	0.0875	0.0273	1.3363	0.0199	0.2405	1.7344	11.5987	0.0121	1.5185
Pine	59.8450	0.0081	0.0860	0.0251	0.8308	0.0174	0.2056	1.3111	11.3134	0.0098	0.9627
Spruce	77.4686	0.0030	0.0836	0.0280	0.6393	0.0250	0.2851	1.7769	26.8149	0.0120	1.3815
Birch	23.6670	0.0234	0.0802	0.0258	1.7425	0.0200	0.1823	2.0643	8.6172	0.0131	0.3988

. All parameters were statistically significant (p < 0.05).

3.2. Bivariate and Marginal Distributions

After the fixed-effect parameters were obtained in Section 3.1, which are listed in Table 2, they were used to evaluate the accuracy of the prediction and forecasting in subsequent sections by using data from the validation dataset. The bivariate mixed-effect parameters SDE model was developed by combining the two univariate models through a bivariate stochastic process. The model considered two correlated observations, tree diameter and polygon area, reflecting the high variation of stand density among stands of Lithuania. The main goal in an SDE modeling framework is to determine the probability density function of the solution, formulated as a diffusion process since it capacitates for calculating any univariate moment, which allows for calculating the mean, variance, and tolerance regions for bivariate cases. In our setting, the newly developed bivariate density and marginal and conditional densities were visually evaluated corresponding to the diameter and polygon area observed data from the validation dataset at a given average stand age. Therefore, the corresponding observed data from the validation dataset and its fitted distribution were graphically visualized. The Voronoi tessellations presented in Figure 1 show dynamic of polygon areas via stand age (time). These diagrams also demonstrate information about variations in the polygon area, which grows with age. Figure 2, Figure 3, Figure 4 and Figure 5 illustrate the fitted marginal probability density functions and underlying frequency distributions (histograms) for new plots from the validation dataset with two cycles of remeasurements at an average stand age. All the fitted probability density functions take the values of the fixed-effect parameters from Table 2. Random effects ${\phi }_d,{\phi }_p,\mathrm{and} {\phi }_0$ for a new plot from the validation dataset were calibrated as:

$$\stackrel{\wedge }{\phi }=\underset{\left({\phi }_d,{\phi }_p,{\phi }_0\right)}{argmax}\left({{\displaystyle \sum }}_{j=1}^{m}ln\left(f\left(x_{1j},x_{2j},t_j\left|\stackrel{\wedge }{\theta },{\phi }_d,{\phi }_p,{\phi }_0\right.\right)\right)+ln\left(\varphi ({\phi }_d\left|{\stackrel{\wedge }{\sigma }}_d)\right.\right)+ln\left(\varphi \left({\phi }_p\left|{\stackrel{\wedge }{\sigma }}_p\right.\right)\right)+ln\left(\varphi \left({\phi }_0\left|{\stackrel{\wedge }{\sigma }}_0\right.\right)\right)\right)$$

(27)

where $\left\{\left(x_{11},x_{21}\right),\left(x_{12},x_{22}\right),\dots ,\left(x_{1m},x_{2m}\right)\right\}$ is the newly observed dataset, t₀ = 4, $\left(x_{10},x_{20}\right)=\left(0.1,\widehat{\delta }+{\phi }_0\right)$, $f\left(x_1,x_2,t\left|\stackrel{\wedge }{{\theta }^1},{\phi }_d,{\phi }_p\right.,{\phi }_0\right)$ is the bivariate normal-lognormal density function with the mean and variance defined by Equations (5) and (6), $\varphi (·\left|\sigma )\right.$ is the univariate normal density function with zero mean and standard deviation σ, and the estimated values of fixed-effect parameters are denoted by “hat” (listed in Table 2).

Figure 2, Figure 3, Figure 4 and Figure 5 show that univariate conditional probability density functions defined by inserting additional explanatory variable polygon areas or diameters undergo indistinguishable changes in comparison with corresponding univariate marginal densities. The similarity of the marginal and conditional distributions suggests that the additional explanatory variable (polygon area or diameter) has a relatively small effect on the values of the response variable (diameter or polygon area). Moreover, the figures superpose frequency distribution and estimated normal or lognormal probability distributions that correspond to the probability density function of the solutions of SDE (2). Frequency distributions are not exactly in agreement with the normal or lognormal distribution shapes due to the variation in tree age or probably did not arise from our presented form. The tree diameter histograms shown in Figure 2, Figure 3, Figure 4 and Figure 5 confirm the assumption that the tree-diameter distribution in a particular forest stand has an approximately symmetrical shape. The lack of symmetry in Figure 4 for Norway spruce trees occurs probably due to the planned thinning in the stand, as lightening works are usually carried out in stands of this age. The symmetry and standard deviation of the diameter distribution increase with age. In contrast, the polygon area histograms reveal asymmetry of the distribution, and the long tail extends to the right. Histograms of the polygon area show a number of unusual observed values. This scenario may be the result of a directional felling or soil properties. Lastly, the newly derived bivariate probability density function was a good match for our validation data.

For the observed diameter and polygon area data fitted with the bivariate Vasicek–Gompertz-type SDE (2), the tolerance region of mean vector ${\mu }^i(t|t_0,x_0)$, defined by Equation (5), takes the following well-known inequality [22]:

$${\left(\left(\begin{array}{c}{x}_1\\ ln\left(x_2\right)\end{array}\right)-{\mu }^i\left(t|t_0,x_0\right)\right)}^T{\left[\Sigma \left(t|t_0\right)\right]}^{-1}\left(\left(\begin{array}{c}{x}_1\\ ln\left(x_2\right)\end{array}\right)-{\mu }^i\left(t|t_0,x_0\right)\right)\le \chi ,$$

(28)

where x is the tolerance coefficient [23]. For lognormally distributed polygon area data, we plot a tolerance region using a logarithmic axis. Figure 6, Figure 7, Figure 8 and Figure 9 show the tolerance regions for β = 0.95 and confidence level γ = 0.95, which correspond to the randomly selected stand from the validation dataset with two cycles of remeasurements. The random effects were calibrated by Equation (27). For the tolerance region plots, the x value was chosen from Table 1 in [23]: the setting β = 0.95 and γ = 0.95 produces a x of 10.02. The bivariate tolerance region for the Vasicek–Gompertz-type model (2) enables us to decide whether the newly developed distribution function corresponds well with the observed data of the diameter and polygon area. Figure 6, Figure 7, Figure 8 and Figure 9 illustrate that the 95% tolerance regions had reasonable coverage rates for different tree-species scenarios. Information concerning tolerance regions is implicit in the hybrid bivariate distribution derived for the Vasicek–Gompertz-type SDE (2).

4. Discussion

Figure 2, Figure 3, Figure 4 and Figure 5 show that very high variation existed in the shape of the diameter and polygon area frequency distributions among plots of a given stand age. In our models, this variation was well accounted for by the used diffusion process defined by SDE (2). Since the observed sample plots represent forest stands in the region, the plot effects were generally included by using three random variables (effects). To properly test our mixed-effect model of diameter, polygon area, and number of trees per hectare, we used observed sample plots from the validation dataset to show its robustness in predicting (via current age) and forecasting (future age).

4.1. Modeling Tree-Diameter Dynamics: Predicting and Forecasting

Traditionally, individual-tree-diameter-growth regression models describe growth as a function of an age (tree or stand). Most individual-tree-diameter-growth models were framed using an algebraic difference approach [24] and its mathematical generalizations [25]. The SDE (2) developed in this study enables us to describe a wide range of tree- and stand-growth variables. To apply the mixed-effect SDE model defined by Equation (2), we had two adaptation strategies—prediction and forecasting. First, the underlying random effects of the model were precisely calibrated (using all remeasurement cycles from the validation dataset). Hence, random effects were calibrated by Equation (27) using data from the validation dataset in order to define the predictions (dynamic against the age) of the mean, variance, and quantiles of the diameter or polygon area in a particular stand. In this strategy, we could average the results of multiple realizations of tree diameters in different plots and obtain important characteristics (mean and variance) that could not be seen in one tree realization. Second, the underlying random effects of the model were not precisely calibrated (using only the first measurement cycle in a plot). Hence, random effects were calibrated by Equation (27) using measurements of diameter, polygon area, and age at the initial measurement cycle from the validation dataset in order to define the forecasts, at the 5-, 13-, and 35-year forecast periods, of the diameter or polygon area for each individual tree. Most remeasured plots from the validation dataset were unevenly aged and mixed-species. Figure and table illustrations are separately presented for all tree species, Scots pine tree species, Norway spruce tree species, and silver birch tree species. Figure 10 shows predictions of the mean diameter, 0.05 and 0.95 quantiles for all, Scots pine, Norway spruce, and silver birch tree species for two randomly selected stands from the validation dataset.

The mean prediction error (percentage of prediction error, %) B, mean absolute prediction error (percentage of absolute prediction error, %) AB, root-mean-square error (percentage of root-mean-square error, %) RMSE, and coefficient of determination R² were used to evaluate the results of the mean diameter-marginal and conditional model fit defined by Equations (7) and (19), respectively. The calculated results of statistical measures using the validation dataset, the fixed-effect parameters from Table 2, and random effects calibrated by Equation (27) are presented in

Table 3. Statistical measures for marginal and conditional models of mean (stand) diameter predictions.

Tree Species	Marginal Mean (Equation (7))				Conditional Mean (Equation (19))
	B (%)	AB (%)	RMSE (%)	R2	B (%)	AB (%)	RMSE (%)	R2
All	0.019 (0.10)	0.798 (3.99)	1.133 (5.67)	0.949	−0.431 (−2.16)	0.782 (3.91)	0.989 (4.95)	0.961
Pine	−0.131 (−0.50)	0.841 (3.19)	1.141 (4.33)	0.973	−0.410 (−1.56)	0.846 (3.21)	1.158 (4.39)	0.972
Spruce	0.130 (0.90)	0.977 (6.80)	1.258 (8.76)	0.928	−0.345 (−2.40)	0.906 (6.30)	1.167 (8.12)	0.938
Birch	−1.038 (−5.94)	1.775 (10.15)	3.096 (17.71)	0.816	−1.354 (−7.74)	1.872 (10.71)	3.086 (17.65)	0.817

. The prediction performance of both models (marginal and conditional) showed that both models were highly capable of identifying the mean value of the diameter in a stand. Subsequently, only a small improvement in the statistical measures was found when we used the polygon area as an additional explanatory variable (conditional model (19)). Table 3 illustrates the accuracy of mean diameter predictions by using the observed validation dataset.

The marginal and conditional models defined by Equations (7) and (19) can also be used successfully as individual-tree-based models in forecasting tree growth regardless of species, age, and tree polygon area. In Equation (7), if initial point $\left(x_{10},x_{20}\right)=\left(0.1,\hat{\delta }+{\phi }_0^i\right)$ was changed by point $\left(x_{10},x_{20}\right)=\left(d_{in,j}^i,p_{in,j}^i+0.\right),$ i = 1,…,K (K is the number of the observed sample plots from the validation dataset), we could calculate forecasts of individual-tree diameters and compare them with the observed dataset by determining the duration of the forecast period (5, 13, and 35 years), where $\left(d_{in,j}^i,p_{in,j}^i\right)$ is the diameter and polygon area of the j-th tree in the i-th plot at base age t_in,j, and random effects ${\phi }_d^i$${\phi }_p^i$ and ${\phi }_0^i$ are equated to 0.

Table 4. Statistical measures of diameter forecasts for 5-, 13-, and 35-year forecast periods.

Tree Species	5-Year Forecast Period				13-Year Forecast Period				35-Year Forecast Period
	B (%)	AB (%)	RMSE (%)	R2	B (%)	AB (%)	RMSE (%)	R2	B (%)	AB (%)	RMSE (%)	R2
All	−0.054 (−0.31)	0.949 (5.52)	1.440 (8.37)	0.977	−0.185 (−094)	2.140 (10.82)	2.911 (14.72)	0.914	−0.359 (−1.43)	4.606 (18.45)	5.654 (22.65)	0.739
Pine	−0.335 (−1.48)	0.982 (4.34)	1.447 (6.39)	0.968	−0.919 (−3.65)	2.178 (8.65)	2.816 (11.19)	0.886	−2.352 (−7.76)	4.469 (14.74)	5.674 (18.71)	0.615
Spruce	0.147 (1.24)	0.880 (7.42)	1.402 (11.82)	0.965	0.255 (1.83)	1.927 (13.85)	3.005 (21.60)	0.859	0.681 (3.62)	4.123 (21.89)	5.159 (27.39)	0.736
Birch	0.020 0.10)	1.286 (6.85)	1.757 (9.38)	0.960	−0.321 (−1.50)	2.697 (12.59)	3.563 (16.64)	0.832	1.289 (4.92)	4.196 (16.03)	5.405 (20.65)	0.679

shows the forecast statistical measures of a tree-individual scenario model calculated for the 5-, 13-, and 35-year forecast periods using the fixed-effect parameters estimates in Table 2.

In general, the relative success of the 35-year forecast period in forecasting individual tree diameters showed that the values of statistical measures were considerably smaller than statistical measures were for the 5- and 13-year forecast periods. Therefore, the shorter (5 and 13 years) forecast periods provided us with reasonably accurate forecasts of tree diameters but performed quite poorly farther into the future. Compared to different tree species, Norway spruce species produced the lowest value of all statistical measures.

4.2. Modeling Tree Polygon Area: Predicting and Forecasting

The tree polygon areas show spatial tessellation based on closeness to trees in a particular stand. It was suggested by mathematician Voronoi [10] and named after him: Voronoi polygons (diagrams). A particular forest stand is driven with occasional events, and the spatial tree arrangement in it is featured with complexity and variability. Most studied Voronoi polygons are only static. A first approach for the dynamization of Voronoi diagrams could be an investigation by a procedure for inserting and deleting single trees (points), each in linear time. In modeling real dynamic scenes of a particular forest stand, parallel continuous changes of trees with a fast update of the Voronoi diagram are desirable. Four realizations at different times of a Voronoi polygon in a particular stand are presented in Figure 1. The main result of the present work consists in the dynamization of the underlying Voronoi polygon areas by the Gompertz-type diffusion process. Figure 11 shows predictions of the mean polygon area and 0.05 and 0.95 quantiles for two randomly selected stands from the validation dataset. Random effects were calibrated by Equation (27) using measurements of tree positions from the validation dataset, and the fixed-effect parameters are from Table 2.

The results of statistical measures of polygon area predictions, calculated using tree positions from the validation dataset, the fixed-effect parameters from Table 2, and random effects calibrated by Equation (27), are presented in

Table 5. Statistical measures for marginal and conditional models of mean polygon area predictions.

Tree Species	Marginal Mean (Equation (10))				Conditional Mean (Equation (22))
	B (%)	AB (%)	RMSE (%)	R2	B (%)	AB (%)	RMSE (%)	R2
All	−0.305 (−2.46)	1.062 (8.58)	1.410 (11.26)	0.929	−0.184 (−1.49)	0.914 (7.39)	1.184 (9.57)	0.949
Pine	−0.392 (−3.17)	0.726 (5.87)	1.009 (8.16)	0.958	−0.276 (−2.23)	0.651 (5.27)	0.883 (7.14)	0.968
Spruce	−0.441 (−3.78)	1.155 (9.89)	1.597 (13.68)	0.913	−0.283 (−2.42)	1.062 (9.10)	1.454 (12.46)	0.928
Birch	0.202 (1.84)	1.628 (14.84)	2.318 (21.13)	0.790	0.442 (4.03)	1.571 (14.32)	2.145 (19.56)	0.820

, which clarifies the importance of the diameter employed in the modeling process. The prediction performance of both models defined by Equations (10) and (22) showed that both models were highly capable of identifying the mean value of the polygon area in a plot. On the other hand, only small improvement in the statistical measures was found (up to 3% in the coefficient of determination) when we used diameter as an additional explanatory variable in the conditional model defined by Equation (22).

Statistical measures of forecasts of trees’ individual polygon areas using an observed validation dataset, the fixed-effect parameters from Table 2, and at random effects equal to zero, ${\phi }_d^i=0., {\phi }_p^i=0., {\phi }_0^i=0.$, i = 1, …, K, are presented in

Table 6. Statistical measures for the marginal model (Equation (10)) of polygon area forecasts for 5-, 13-, and 35-year forecast periods.

Tree Species	5-Year Forecast Period				13-Year Forecast Period				35-Year Forecast Period
	B (%)	AB (%)	RMSE (%)	R2	B (%)	AB (%)	RMSE (%)	R2	B (%)	AB (%)	RMSE (%)	R2
All	−0.601 (−6.01)	1.126 (11.26)	1.539 (15.38)	0.959	−0.555 (−4.69)	2.756 (23.31)	3.963 (33.51)	0.786	−1.581 (−11.09)	5.107 (35.83)	6.746 (47.33)	0.463
Pine	−0.762 (−6.34)	1.253 (10.43)	1.621 (13.49)	0.959	−1.194 (−8.72)	2.969 (21.70)	3.850 (28.13)	0.803	−3.655 (−23.78)	5.959 (38.66)	7.167 (46.50)	0.408
Spruce	−0.384 (−4.65)	0.983 (11.90)	1.466 (17.74)	0.953	0.221 (2.19)	2.601 (25.81)	4.238 (42.04)	0.725	0.664 (5.12)	4.575 (35.25)	6.798 (52.38)	0.405
Birch	−0.539 (−6.04)	1.004 (11.25)	1.198 (13.43)	0.971	−0.785 (−7.93)	2.504 (25.31)	3.373 (34.09)	0.827	0.201 (1.37)	5.788 (39.61)	7.458 (51.13)	0.491

for the forecast periods of 5, 13, and 35 years. The conditional model defined by Equation (22) was evaluated using a data subset of diameters drawn from the validation dataset at the projected (forecast period) age and showed that the influence of diameters cannot greatly improve the forecast ability of a mixed-effect polygon area model. Statistical measures in Table 6 show that the accuracy of the model forecast decreases significantly with increasing the forecast period to 35 years.

4.3. Modeling Stand Density: Predicting and Forecasting

To manage the evolution of the number of trees per hectare from the early sapling stage to any stage in mixed-species, unevenly aged forests, reliable predictive and forecast models are needed. The complete size–density trajectory of a stand from an early development stage follows the form framed by the maximal size–density relationship [26]. Recently, diffusion processes were used to define maximal size–density equations [27]. Many various stand-density measures were developed as relationships of mean area available to trees in a particular stand [28,29]. The stand density expresses a stand occupancy in abstract form; consequently, in this study, the stand density per hectare dynamic is related to the dynamic of the polygon area, defined by Equations (10) and (22), in the following forms for all the tree species:

$$N^i\left(t|t_0,x_{20}\right)=\frac{10000}{m_2^i\left(t|t_0,x_{20}\right)},$$

(29)

$$N^i\left(t|t_0,x_{20},x_{11}\right)=\frac{10000}{{\eta }_2^i\left(t,x_{11}|t_0,x_0\right)},$$

(30)

and for constituent tree species:

$$N^i\left(t|t_0,x_{20}\right)=o{c}_{in}\frac{10000}{m_2^i\left(t|t_0,x_{20}\right)},$$

(31)

$$N^i\left(t|t_0,x_{20},x_{11}\right)=o{c}_{in}\frac{10000}{{\eta }_2^i\left(t,x_{11}|t_0,x_0\right)},$$

(32)

where $o{c}_{in}$ is the occupation proportion of a specific tree species in a stand at an age of the first measurement (0 < $o{c}_{in}$ < 1).

Sustainable forest management requires the comprehensive understanding of the long-term dynamics of stand density. The stand-density models defined by Equations (29)–(32) enable us to evaluate stand density in both prediction and forecasting scenarios. For calculating predictions of stand density, we used the fixed-effect parameters from Table 2 and the random effects calibrated by Equation (27), using a full validation dataset. Figure 12 shows the mean stand-density dynamics for all tree species and constituent tree species scenarios and compares with observed datasets for three randomly selected stands from the validation dataset. The accuracy measures of the mean stand-density predictions are presented in

Table 7. Statistical measures for marginal and conditional models of mean stand-density predictions.

Tree Species	Marginal Mean (Equation (29) or (31))				Conditional Mean (Equation (30) or (32))
	B (%)	AB (%)	RMSE (%)	R2	B (%)	AB (%)	RMSE (%)	R2
All	0.686 (0.07)	102.23 (10.50)	166.149 (17.07)	0.873	−8.253 (−0.84)	90.608 (9.31)	147.297 (15.13)	0.901
Pine	19.339 (4.30)	37.682 (8.39)	49.302 (10.97)	0.966	16.429 (3.65)	34.762 (7.74)	44.767 (9.96)	0.972
Spruce	−3.720 (−0.62)	76.165 (12.74)	123.67 (20.67)	0.934	−14.901 (−2.49)	68.768 (11.50)	110.925 (18.56)	0.947
Birch	5.177 (7.89)	9.920 (15.12)	14.118 (21.52)	0.984	4.860 (7.41)	10.214 (15.57)	14.151 (21.57)	0.984

.

The Voronoi polygon of tree positions is a powerful tool for understanding the spatial competition properties on the basis of closeness to trees in a particular stand. Relatively few studies reported the construction of Voronoi polygons in forest-stand modeling [30] despite a wide array of potential applications. Dynamic Voronoi polygons have applications among stand-density models, as they can be used to accurately define the rate of natural mortality among trees within a forest. Results of diameter importance for stand-density modeling processes are shown in the conditional models defined by Equations (30) and (32), presented in Table 7, which shows that employing the diameter in the modeling process provided small improvement in statistical measures (up to 3% in the coefficient of determination).

In stand-density dynamic models, transition probability density functions only consider natural mortality that is influenced by competition, whereas general response functions account for the effects of the site index, basal area, and other factors [31]. The newly developed stand-density dynamic functions defined by Equations (29)–(32) were developed with the purpose of formalizing a whole-stand model for predicting and forecasting a given mixed-species forest stand. As such, the observed dataset (40 plots) from long-term remeasurements (from 1 until 6 cycles) of permanent plots was used to evaluate the fixed-effect parameters for both predicting and forecasting scenarios. Concerning the stand-density forecast scenario, the random effects for a new stand were calibrated by Equation (27) using the first measurement cycle observed sample from the validation dataset. Comparisons of stand-density forecasts (projections) among the forecast periods of 5, 13, and 35 years are presented in

Table 8. Statistical measures of mean stand-density forecasts (Equations (29) and (31)) for 5-, 13-, and 35-year forecast periods.

Tree Species	5-Year Forecast Period				13-Year Forecast Period				35-Year Forecast Period
	B (%)	AB (%)	RMSE (%)	R2	B (%)	AB (%)	RMSE (%)	R2	B (%)	AB (%)	RMSE (%)	R2
All	−59.492 (−5.36)	60.113 (55.41)	72.799 (6.56)	0.981	−24.644 (−2.73)	94.991 (10.54)	118.212 (13.11)	0.897	62.025 (8.53)	103.635 (14.67)	125.749 (17.30)	0.605
Pine	−15.920 (−3.25)	40.292 (8.24)	53.640 (10.97)	0.966	31.761 (7.36)	46.707 (10.82)	56.766 (13.15)	0.938	74.438 (20.26)	87.838 (23.90)	105.201 (28.63)	0.751
Spruce	−24.434 (−5.50)	37.735 (8.50)	56.291 (12.68)	0.982	−58.010 (−15.04)	60.760 (15.73)	104.466 (27.08)	0.914	−19.460 (−6.57)	29.808 (23.90)	51.048 (17.26)	0.925
Birch	−68.224 (−9.68)	70.412 (9.99)	90.988 (12.91)	0.973	−82.303 (−15.26)	90.418 (16.76)	128.911 (23.90)	0.891	−36.060 (−8.68)	59.143 (14.25)	78.997 (19.03)	0.900

.

All the tested models for all the tree species provided good forecasts for data with high statistical measures. All the models for Scots pine, Norway spruce, silver birch, and all the tree species resulted in high coefficients of variation for the forecast periods of 5, 13, and 35 years in intervals of 96.6–98.2%, 89.1–93.8%, and 60.5–92.5%, respectively. Both scenario models for Scots pine trees were comparably better. Similar results were obtained in model ranking on the basis of B, B%, AB, AB%, RMSE, and RMSE%. Statistically, the marginal models defined by Equations (29) and (30) and the conditional models defined by Equations (30) and (32) provided similar forecasts for all forecast periods.

5. Conclusions

In this paper, we have studied the evolution of the mixed-species, unevenly aged forest stands by using a bivariate hybrid diffusion process and Voronoi diagram. The growth model considers two different system states (tree diameter and polygon area). In summary, derived individual-tree and whole-stand models describe how trees grow in diameter and how forest-stand structures are modified over time. On the other hand, one of the most fundamental features of our developed models is that they are strongly symmetrical, as they allow for forecasting trajectories in the future and the past. Numerical example by using experimental sample plots in Lithuania with measurements of tree position, age, and diameter at breast height showed high accuracy of the obtained results and the importance of the work. From a statistical point of view, the newly developed growth models produced acceptable predictions and forecasts. Our proposed bivariate hybrid diffusion process and Voronoi diagram approach also outperformed other existing techniques [32,33]. The newly developed hybrid bivariate distribution also provides a further in-depth understanding of the behavior of the stand basal area and volume.

Future work should try to extend our work to describe hybrid 3-, 4-, and 5-variate diffusion processes and copula approach for developing the link between state variables, as an example, diameter, height, crown width, crown base height, and available polygon area.