The horizontal structure of a forest stand can be characterized by the distribution of tree diameter at breast height (at 1.37 m, DBH) and the relative locations of trees, while the vertical structure of the same stand can be represented by the distribution of tree height (HT) and diameter-height relationships. Various probability density functions (pdf) have been used to describe the frequency distributions of tree diameter, such as Weibull (Bailey et al., 1973), lognormal (Bliss et al., 1964), gamma (Nelson, 1964), beta (Clutter et al., 1965), Johnson's S_B (Hafley et al., 1977), and logit-logistic (Wang et al., 2005). However, few efforts have been made to fit the frequency distributions of tree height (Schreuder et al., 1977). Thus, in modeling the vertical structure of a stand, the most common practice is to fit a distribution function to the diameter frequency data, and then use an empirical diameter-height relationship to estimate the mean height per diameter class. Although this approach is satisfactory in many situations, it ignores the natural relationships between tree diameter and height by treating them separately since the tree height can vary considerably for a given diameter (Scheuder et al., 1977; Tewari et al., 1999).

An alternative is to model the distributions of tree diameter and height simultaneously using a bivariate distribution. This approach is desirable because it allows us to predict the tree height by diameter from the fitted bivariate distribution through the conditional expectation of tree height given DBH. Over the last decades, Johnson's S_BB have been utilized for modeling the bivariate tree diameter-height frequency data (Schreuder et al., 1977; Hafley et al., 1985; Knoebel et al., 1991; Tewari et al., 1999; Zucchini et al., 2001). Li et al. (2002) applied a bivariate generalized beta distribution (GBD-2) to model the joint distributions of tree diameter and height of Douglas-fir stands in the northwest, USA, and found that GBD-2 performed better than Johnson's S_BB in goodness-of-fit statistics and stand volume predictions. Wang et al.(2007; 2008; 2010) used a general copula approach to modeling multivariate distributions, with numerical examples for bivariate tree diameter-height data or trivariate tree diameter-height-volume data.

The purposes of this study were 1) to fit two bivariate distribution functions (i. e., GBD-2 and Johnson's S_BB) to the joint distributions of tree diameter and height of spruce-fir stands in the northeast, USA, 2) to compare the model fitting and performance between the two bivariate distribution functions, and 3) more importantly, to develop regression models for predicting the parameters of the GBD-2 distribution using available stand variables.

Fifty (50) plots were collected from the evenaged, unmanaged natural spruce-fir forests in northwestern Maine, USA, located in the region between 69°W and 71°W in longitude, between 45°N and 46.5°N in latitude and between 750 m and 1 200 m in elevation (Kleinschmidt et al., 1980). The plot area ranged from 0.002 5 to 0.02 hm² in size (mean plot size was 0.008 hm² with stand ard deviation 0.004 hm²). In these plots, balsam fir (Abies balsamea) and red spruce (Picea rubens) accounted for about 95% of total number of trees and 94% of total volume. Other minor species included black cherry (Prunus serotina), eastern white pine (Pinus strobus), white spruce (Picea glauca), black spruce (Picea mariana), and other hardwoods. In each plot, tree DBH, H, crown length (top to the base of crown), and average crown width were recorded for each tree (Solomon et al., 2002). The numbers of trees in each plot ranged from 12 to 109. Mean tree diameters were from 2.2 to 19.2 cm and mean total heights were from 2.7 to 17.1 m (Tab. 1).

Tab.1 Summary statistics of the stand variables across the 50 plots

Tab.1 Summary statistics of the stand variables across the 50 plots Stand variables Mean Std. Min. Max. Number of trees per plot 34.4 17.7 12 109 Stand density/(tree·m-2) 0.61 0.54 0.18 2.52 Stand basal area/(m2·hm-2) 50.2 14.3 10.0 79.6 Stand mean DBH/cm 11.5 3.9 2.2 19.2 Stand mean height/m 11.6 3.7 2.7 17.1 Stand mean crown length/m 4.4 1.5 1.6 8.4 Stand mean crown width/m 1.1 0.3 0.5 1.9

The plot of "Solution Available" region was constructed for the 50 spruce-fir plots (Li et al., 2002). Most of the tree diameter and height distributions of all 50 plots fell in the region of generalized beta distribution (GBD) or the overlap region of GBD and generalized lambda distribution (GLD). Therefore, GBD is appropriate to fit the marginal distributions of tree diameter and height for the 50 plots and GBD-2 is applied to model the joint distribution of tree diameter and height.

The functions and methods in Li et al. (2002) were used to 1) fit the marginal GBD distributions to tree diameter and height, respectively; 2) fit GBD-2 to the joint distribution of tree diameter and height. The probability density functions of GBD and GBD-2 are provided in Appendix 1; 3) fit the marginal S_B distributions to tree diameter and height, respectively; 4) fit Johnson's S_BB to the joint distribution of tree diameter and height. The probability density functions of S_B and S_BB are provided in Appendix 2; 5) test the goodness-of-fit for the marginal GBD or Johnson's S_B distributions of tree diameter and height, as well as the two bivariate distributions; and 6) compare tree height and volume predictions between GBD-2 and Johnson's S_BB (Appendix 4). Prediction bias was defined as the difference between observed and predicted tree height or volumes, and relative bias was calculated as the percentage of the prediction bias over the predicted tree height or volume, respectively. An average relative bias was then obtained for each of the 50 plots. To examine the prediction bias across tree sizes, all trees in the 50 plots were grouped into 2-cm diameter classes. An average relative bias was calculated for each diameter class.

In forestry practice, it is desirable to be able to predict the parameters of the GBD-2 distributions based on conventional forest inventory variables. One intuitive approach is to build empirical regression models between these parameters and stand variables. However, the correlations between the parameters of a frequency distribution and stand variables are generally low due to the fact that parameter estimates are highly data-specific, stands with similar conditions can produce highly variable parameter estimates. Therefore, it is difficult to directly predict the parameters based on stand characteristics (Schreuder et al., 1977; Knoebel et al., 1991). Alternatively, parameter recovery approach can be employed, i. e., instead of predicting the GBD-2 parameters directly from stand variables, some moments of tree diameter and height can be predicted from the stand variables because these moments are highly correlated with other stand variables. In this study, the seven linear regression models were developed to predict the response variable Y which included ψ (an association measure between two marginal distributions), minimum of DBH, maximum of DBH, stand ard deviation of DBH, minimum of H, maximum of H, and stand ard deviation of H:

where α₀-α₆ are regression coefficients to be estimated, and ε is the model error term. Consequently, the parameters of GBD-2 can then be solved by the relationships among the parameters using equations (8) - (11) in Appendix 3 (Li et al., 2002).

Of the 50 plots none (0%) of the diameter distributions and one (2.0%) of the height distributions fitted by the GBD distribution were significantly different from the observed distribution at α = 0.05 level. Correspondingly, three (6.0%) of the diameter distributions and seven (14.0%) of the height distributions fitted by the Johnson's S_B distribution were significantly different from the observed distribution. Further, none (0%) of the predicted bivariate GBD-2 distributions was significantly different from the observed distributions at α = 0.05 level, while nine (18.0%) of the predicted Johnson's S_BB distributions were significantly different from the observed ones. Therefore, the GBD-2 distribution was better than the Johnson's S_BB distribution in fitting the joint distribution of tree diameter and height for the spruce-fir stands in the Northeast, USA.

For predicting tree height from the two joint distributions, the predicted tree height using GBD-2 was significantly different (α = 0.05) from the observed tree height in only 1 (2.0%) of the 50 plots, but this occurred in 9 (18.0%) of the 50 plots for the Johnson's S_BB . Again, the GBD-2 performed better in predicting tree height, although numerical methods must be employed for its use.

Fig. 1 shows the average relative biases (%) of volume predictions across 2 cm diameter classes. It appears that the Johnson's S_BB distribution produces relatively larger positive biases (underestimation) for small-sized trees (e. g., < 10 cm in diameter) and for large-sized trees (e. g., > 22 cm in diameter), but relatively larger negative biases (overestimation) for middle-sized trees than did the GBD-2 distribution.

	Fig. 1 Average relative bias (%) of volume prediction from the GBD-2 and Johnson's SBB distributions for diameter classes over all of the 50 plots

The correlations between the GBD-2 parameters and available stand variables in the data were low in general, and about 52% (28 /54) of these correlation coefficients were not significantly different from zero at α = 0.05 level (Tab. 2). On the other h and, the estimated association parameter $\mathop \psi \limits^ \wedge $ of GBD-2, the minimum values, maximum values and the stand ard deviation of tree diameter and height were significantly correlated with majority (88% or 37/42) of the available stand variables (Tab. 3). Therefore, seven linear regression models for the $\mathop \psi \limits^ \wedge $, the minimum values, maximum values and the stand ard deviation of tree diameter and height were constructed using these stand variables including stand density, basal area, mean tree diameter, mean tree height, and mean crown length and width. Since the stand variables themselves were highly correlated, the principal component regression approach was employed (Neter et al., 1990). Tab. 4 provides the final regression coefficients and model fitting statistics for the seven regression models. The model R² were ranged from 0.5 to 0.95 for the models developed for predicting the minimum and maximum values and the stand ard deviation of tree diameter and height, while the model R² for the model of the parameter $\mathop \psi \limits^ \wedge $ of the GBD-2 was relatively low (0.28).

Tab.2 Correlations of the parameters of GBD-2 with the available stand variables

Tab.2 Correlations of the parameters of GBD-2 with the available stand variables GBD-2 parameters Density /(tree·m-2) Basal area /(m2·hm-2) Mean DBH /cm Mean height /m Mean crown length /m Mean crown width /m $\mathop \psi \limits^ \wedge $(P-value) 0.344 9(0.012 3) -0.105 8(0.455 3) -0.413 0(0.002 3) -0.398 6(0.003 4) -0.392 8(0.004 0) -0.167 3(0.255 7) DBH ${\mathop \beta \limits^ \wedge _1}$(P-value) -0.339 9(0.013 7) 0.163 2(0.247 7) 0.354 4(0.009 9) 0.421 6(0.001 9) 0.299 8(0.030 8) 0.131 0(0.374 8) ${\mathop \beta \limits^ \wedge _2}$(P-value) -0.414 4(0.002 3) 0.544 4(0.000 0) 0.462 7(0.000 6) 0.391 1(0.004 1) 0.398 4(0.003 4) 0.061 0(0.680 4) ${\mathop \beta \limits^ \wedge _3}$(P-value) -0.204 9(0.145 1) 0.154 8(0.273 2) 0.248 3(0.075 9) 0.224 8(0.109 1) 0.067 3(0.635 4) 0.047 4(0.749 1) ${\mathop \beta \limits^ \wedge _4}$(P-value) 0.031 5(0.824 6) 0.048 1(0.734 9) -0.049 4(0.728 1) -0.050 6(0.721 6) -0.088 3(0.533 5) -0.258 8(0.075 7) Height ${\mathop \beta \limits^ \wedge _1}$(P-value) -0.156 5(0.267 8) 0.093 5(0.509 6) 0.095 4(0.501 2) 0.229 6(0.101 6) 0.127 7(0.366 9) -0.048 6(0.742 7) ${\mathop \beta \limits^ \wedge _2}$(P-value) -0.442 5(0.001 0) 0.364 4(0.007 9) 0.526 4(0.000 1) 0.416 3(0.002 1) 0.332 4(0.016 1) 0.335 0(0.020 0) ${\mathop \beta \limits^ \wedge _3}$(P-value) -0.261 0(0.061 6) 0.099 3(0.483 7) 0.344 0(0.012 5) 0.311 0(0.024 8) 0.077 3(0.586 1) 0.183 7(0.211 2) ${\mathop \beta \limits^ \wedge _4}$(P-value) 0.379 6(0.005 5) -0.378 8(0.005 6) -0.420 5(0.001 9) -0.350 7(0.010 8) -0.464 0(0.000 5) -0.092 6(0.531 4)

Tab.3 Correlations of the parameter $\mathop \psi \limits^ \wedge $ of GBD-2，the minimum values，maximum values and the standard deviations of tree diameter and heights with the available stand variables

Tab.3 Correlations of the parameter $\mathop \psi \limits^ \wedge $ of GBD-2，the minimum values，maximum values and the standard deviations of tree diameter and heights with the available stand variables GBD-2 parameters Density /(tree·m-2) Basal area /(m2·hm-2) Mean DBH /cm Mean height /m Mean crown length /m Mean crown width /m $\mathop \psi \limits^ \wedge $(P-value) 0.344 9(0.012 3) -0.105 8(0.455 3) -0.413 0(0.002 3) -0.398 6(0.003 4) -0.392 8(0.004 0) -0.167 3(0.255 7) DBH Minimum(P-value) -0.682 7(0.000 0) 0.497 2(0.000 2) 0.797 2(0.000 0) 0.838 5(0.000 0) 0.549 1(0.000 0) 0.125 1(0.396 9) Maximum(P-value) -0.782 4(0.000 0) 0.770 4(0.000 0) 0.887 7(0.000 0) 0.832 8(0.000 0) 0.695 6(0.000 0) 0.392 0(0.005 9) StDev(P-value) -0.639 0(0.000 0) 0.727 9(0.000 0) 0.709 8(0.000 0) 0.622 7(0.000 0) 0.651 5(0.000 0) 0.349 5(0.014 9) Height Minimum(P-value) -0.537 6(0.000 0) 0.338 6(0.014 1) 0.632 5(0.000 0) 0.719 6(0.000 0) 0.358 5(0.009 1) 0.037 1(0.802 3) Maximum(P-value) -0.917 7(0.000 0) 0.678 8(0.000 0) 0.933 2(0.000 0) 0.955 6(0.000 0) 0.689 0(0.000 0) 0.427 8(0.002 4) StDev(P-value) -0.445 7(0.000 9) 0.512 0(0.000 1) 0.395 9(0.003 7) 0.318 5(0.021 4) 0.448 9(0.000 8) 0.278 7(0.055 1)

Tab.4 Regression coefficients for the parameter $\mathop \psi \limits^ \wedge $ of GBD-2，the minimum (Min.) values，maximum (Max.) values and the standard deviation (Std.) of tree diameter (DBH) and height (H) using the available stand variables

Tab.4 Regression coefficients for the parameter $\mathop \psi \limits^ \wedge $ of GBD-2，the minimum (Min.) values，maximum (Max.) values and the standard deviation (Std.) of tree diameter (DBH) and height (H) using the available stand variables Dependent variable Intercept Density/(tree·m-2) Basal area /(m2·hm-2) Mean DBH /cm Mean H/m Mean crown length/m Mean crown width/m Fit statistics R2 P-value $\mathop \psi \limits^ \wedge $ 128.963 5 -18.546 8 1.363 4 -3.836 9 -4.067 8 -9.662 4 -2.543 2 0.277 0 0.030 6 Min. H -0.621 5 1.369 7 -0.059 5 0.524 1 0.637 9 -0.488 1 -2.970 9 0.695 9 < 0.000 1 Max. H 8.676 6 -2.268 5 0.023 0 0.279 9 0.344 4 -0.033 9 0.213 2 0.951 1 < 0.000 1 Std. H 2.659 2 -1.085 5 0.034 5 -0.085 4 -0.102 1 0.144 8 0.329 2 0.504 2 < 0.000 1 Min. DBH -1.062 0 0.535 0 -0.017 3 0.329 7 0.388 0 -0.069 1 -1.825 4 0.769 3 < 0.000 1 Max. DBH 2.638 0 -1.565 1 0.159 7 0.345 8 0.345 1 0.324 6 1.549 1 0.868 4 < 0.000 1 Std. DBH 1.117 7 -0.581 2 0.046 0 -0.006 0 -0.024 0 0.183 9 0.463 0 0.750 3 < 0.000 1

To evaluate the performance of the prediction models, the univariate and bivariate χ² tests (Li et al., 2002) were used for testing the goodness-of-fit between the observed and the predicted marginal and joint distributions of tree diameter and height of each plot. The results indicated that 69% of the plots for the diameter distribution, 63% of the plots for the height distribution, and 67% of the plots for the bivariate distribution were predicted satisfactorily by the prediction models.

Three plots were arbitrarily selected as the examples to show the application of the prediction models. The estimated and predicted model parameters as well as the stand variables of these three plots were listed in Tab. 5. The estimated model parameters were obtained from fitting the GBD-2 distribution to each plot, while the predicted model parameters were obtained from the seven prediction models using the available values of the six predictor stand variables of each plot. The scatter plots of the observed and the simulated tree height versus diameter, as well as the histograms of the observed and the simulated tree diameter and height are shown in Fig. 2, 3 and 4. It appears that for all three example plots the predicted height-diameter relationships cover the similar ranges of the observed ones, and the predicted diameter and height frequency histograms have similar patterns and trends to the observed ones.

Tab.5 Estimated (Est.) and predicted (Pred.) parameters of the regression model and the stand variables for the three example plots

Tab.5 Estimated (Est.) and predicted (Pred.) parameters of the regression model and the stand variables for the three example plots Plot 5 26 57 Stand variables Density 0.720 0 0.402 4 0.755 1 Basal area 46.895 2 62.072 8 57.645 7 Mean DBH 8.553 5 13.035 0 9.140 5 Mean H 9.607 4 12.911 0 9.526 8 Mean crown length 3.128 8 5.268 5 2.930 8 Mean crown width 1.114 4 1.389 0 1.194 4 GBD-2 Parameter Est. Pred. Est. Pred. Est. Pred. $\mathop \psi \limits^ \wedge $ 64.000 0 74.581 8 81.000 0 49.158 1 68.000 0 88.372 3 ${\mathop \beta \limits^ \wedge _2}$-DBH 3.250 0 2.809 4 5.364 3 4.487 4 2.900 0 2.672 3 ${\mathop \beta \limits^ \wedge _2}$-DBH 13.700 0 15.204 6 25.027 6 20.256 8 15.000 0 17.238 0 ${\mathop \beta \limits^ \wedge _3}$-DBH 0.377 8 0.149 2 0.234 1 0.291 8 0.249 9 0.149 8 ${\mathop \beta \limits^ \wedge _3}$-DBH 1.181 4 0.892 7 1.792 5 0.769 6 0.754 3 0.914 4 ${\mathop \beta \limits^ \wedge _1}$-H 3.350 0 3.348 3 6.370 0 4.606 5 3.700 0 2.872 2 ${\mathop \beta \limits^ \wedge _2}$-H 10.140 0 10.609 8 12.620 0 12.799 9 10.320 0 11.414 2 ${\mathop \beta \limits^ \wedge _3}$-H 1.368 3 0.830 7 1.929 2 0.772 7 0.462 5 0.643 3 ${\mathop \beta \limits^ \wedge _4}$-H 0.469 5 0.272 6 1.722 3 － 0.040 4 0.127 8 0.175 4

	Fig. 2 Observed versus predicted tree diameter and height of plot 05

	Fig. 3 Observed versus predicted tree diameter and height of plot 26

	Fig. 4 Observed versus predicted tree diameter and height of plot 57

The bivariate distributions of generalized beta distribution (GBD-2) and Johnson's S_BB were fit to a wide range of tree diameter and height distributions of the even-aged, unmanaged natural forest stands of spruce-fir in the Northeast, USA. The goodness-of-fit tests indicated that GBD-2 performed better than did Johnson's S_BB in fitting both marginal and joint distributions of tree diameter and height, as well as in predicting tree height and volume. However, the better predictions of tree height and volume by GBD-2 over S_BB may be partially due to the nature of the regression relationships, i. e., the expectation-type for GBD-2 and median-type for S_BB . Nevertheless, the results in this study were very similar to the previous study for the Douglas-fir stands in the northwest, USA (Li et al., 2002). The weakness of the using the GBD-2 distribution is that the relationship between tree diameter and height is less obvious than the median regression relationship constructed by the Johnson's S_BB distribution. In addition, numerical methods are required for fitting the frequency distributions and predicting tree height using the GBD-2 distribution.

Other shortcomings of this study were 1) the size of the available spruce-fir plots was relatively small. Thus, the parameters of the bivariate distributions may be estimated poorly for the plots with small number of trees, and 2) there were no independent data or plots available for validating the empirical regression models for predicting the model parameters. Therefore, more and large forest plots may be needed to further confirm and validate the results in this study.

In this study empirical regression models were developed for predicting the parameters of the GBD-2 distributions using ordinary stand variables as predictors, such as stand density, basal area, mean tree diameter and total height, and mean crown length and width. Thus, the future stand vertical structure can be predicted when the future values of these stand variables are available.

The major advantage of the bivariate distributions is that it allows for a consistent generation of a bivariate diameter-height system, which can be used to improve the prediction of tree heights and consequently tree volumes. In recent years spatial models have been developed to describe the spatial distribution patterns of trees and stand canopy structure. When modeling three-dimensional profile and changes over time of a forest stand, a bivariate distribution of tree diameter and height is desirable because it can provide a realistic vertical profile and diameter-height dynamics over time.

The probability density function (pdf) of the GBD r and om variable x (e.g., tree diameter or height in this study), if for the parameters β₁, β₂, β₃, β₄ >-1, is

on the interval (β₁, β₁ + β₂), and 0 otherwise, where $C = \frac{{\Gamma \left ({{\beta _3} + {\beta _4} + 2} \right) }}{{\Gamma \left ({{\beta _3} + 1} \right) \Gamma \left ({{\beta _4} + 1} \right) }}$.

The distribution function (H) of the bivariate generalized beta distribution (GBD-2) does not have a closed form, and can be constructed by Plackett's copula (Plackett, 1965) as follows:

where F = F ( x) and G = G (y) are the GBD marginal cumulative distribution functions (cdf) for the two r and om variables, x (e.g., tree diameter) and y (e.g., tree height), respectively, S = 1 + (F + G) (ψ-1), and ψ∈[0, ∞) is a measure of association between the two marginal distributions.Thus, the pdf of GBD-2 is

where f = f (x) and g = g (y) are the GBD marginal pdfs of the two r and om variables, respectively (Li et al., 2002).

The pdf of a univariate Johnson's S_B distribution (Johnson, 1949a) for a r and om variable x is given by

where ε < x < ε + λ, δ > 0, - ∞ < γ, ε < ∞, λ > 0, and ${Z_x} = \gamma + \delta \ln \left ({\frac{{x - \varepsilon }}{{\varepsilon + \lambda - x}}} \right) - N\left ({0, 1} \right) $ 0, 1), with ε the minimum value of x, λ the range of x, γ and δ the shape parameters. From the definition, we know that Johnson's S_B distribution is obtained from the transformation on a st and ard normal variate.

and

where Z₁ and Z₂ have the joint bivariate normal distribution with correlation ρ, and the joint pdf is

One way for estimating the four parameters (β1, β₂, β₃, β₄) of the marginal GBD distribution is as follows: set

and

where x_min and x_max are the minimum values and maximum values of the r and om variable x (e. g., tree diameter or height in this study), respectively. Then, β₃ and β₄ can be estimated as follows:

where $\overline x $ and s² are the sample mean and variance of the GBD r and om variable x, respectively (Li et al., 2002).

For the Johnson's S_BB distribution, the two marginal variables, diameter (x) and height (y), have a simple median regression, namely,

where $\theta = \exp \left[ {\mathop \rho \limits^ \wedge \mathop \gamma \limits^ \wedge } \right]$ and $\phi = \frac{{\mathop \rho \limits^ \wedge {{\mathop \delta \limits^ \wedge }_x}}}{{{{\mathop \delta \limits^ \wedge }_y}}}, {\mathop \varepsilon \limits^ \wedge _x}, {\mathop \lambda \limits^ \wedge _x}, {\mathop \delta \limits^ \wedge _x}, {\mathop \gamma \limits^ \wedge _x}\;\;{\rm{ and }}\;{\mathop \varepsilon \limits^ \wedge _y}, {\mathop \lambda \limits^ \wedge _y}, {\mathop \delta \limits^ \wedge _y}, {\mathop \gamma \limits^ \wedge _y}$ are the estimated parameters of the tree diameter and height distributions, respectively. Hence, this relationship can be naturally used to predict the tree height for a given diameter.

For the GBD-2 distribution, there is no simple relationship between the two marginal variables available. After the parameters are estimated, the conditional expectation of a tree height can be calculated for a given diameter using numerical methods. The expectation of tree height (y) given diameter (x) is

${\rm{where}}\;{g_Y}\left (y \right) = C\beta _{2y}^{ - \left ({{\beta _{3y}} + {\beta _{4y}} + 1} \right) }{\left ({y - {\beta _{1y}}} \right) ^{{\beta _{3y}}}}{\left ({{\beta _{1y}} + {\beta _{2y}} - y} \right) ^{{\beta _{4y}}}}$ is the marginal pdf of the tree height ( y), with $C = \frac{{\Gamma \left ({{\beta _{1y}} + {\beta _{2y}} + 2} \right) }}{{\Gamma \left ({{\beta _{1y}}} \right) \Gamma \left ({{\beta _{2y}}} \right) }}$, S=1+ (F +G) (ψ+1), and F and G are the cumulative density functions of the two marginal variables x and y (Li et al., 2002).