United States

Department of

Agriculture

Forest Service

Southern

Research Station

Research Paper

SRS–31

Recovering Parameters of Johnson’s

SB Distribution

Bernard R. Parresol
Jennifer D. Knoepp, Research Soil Scientist, and

James M. Vose, Project Leader, U.S. Department of

Agriculture, Forest Service, Southern Research Station,

Coweeta Hydrologic Laboratory, 3160 Coweeta Lab Road,

Otto, NC 28763.

The Author

Bernard R. Parresol, Mathematical Statistician, Southern

Forest Inventory, Monitoring, and Analysis Program, U.S.

Department of Agriculture, Forest Service, Southern Research

Station, Asheville, NC 28804.

April 2003

Southern Research Station

P.O.Box 2680

Asheville, NC 28802
1

Recovering Parameters of Johnson’s SB Distribution

Bernard R. Parresol

Abstract

A new parameter recovery model for Johnson’s SB distribution is

developed. This latest alternative approach permits recovery of the range

and both shape parameters. Previous models recovered only the two

shape parameters. Also, a simple procedure for estimating the

distribution minimum from sample values is presented. The new

methodology employs the median and the first and second moments of

the distribution. The methodology is demonstrated by modeling diameter

distributions of unthinned loblolly pine plantations. Compatible

equations for projecting per-hectare values of number of trees and basal

area from initial stand conditions are presented, as well as equations for

predicting median diameter, mean diameter, and the location parameter.

Given estimates of these five stand attributes, the range and the two

shape parameters of the SB distribution can be recovered. The 2 χ

goodness-of-fit test rejected 56 cases out of 527 for conformance to an SB

distribution. Though the SB distribution is very flexible in terms of

distribution shape, about 10 percent of the loblolly plantation

observations did not follow this distribution. Nonetheless, deviation

analysis showed reasonable results, with 77 percent of the variation

explained in current and projected distributions (numbers of trees by

2.5-cm diameter class). Overall, the recovered SB distributions provided

good approximations of the observed diameter distributions.

Keywords: Diameter distributions, kurtosis, Newton-Raphson method,

moments, parameter recovery model, percentiles, Pinus taeda, skewness.

Introduction

Forest managers have long relied on growth-and-yield

forecasts in making management decisions. Growth-andyield models that provide detailed stand distributional

information are particularly useful. Clutter and Bennett

(1965) introduced diameter distribution methodology using

the four-parameter beta probability density function (PDF)

to describe the distribution of the number of trees per unit

area by diameter at breast height (d.b.h.) class. Gove and

Patil (1998) refer to this as the d.b.h.-frequency distribution.

Bailey and Dell (1973) were the first to use the Weibull

probability density function as a d.b.h.-frequency

distribution model. They justified use of the Weibull over

the beta distribution because the Weibull has a closed form

expression for its cumulative distribution function (CDF)

and only three parameters. Hafley and Schreuder (1977)

introduced the SB (system bounded) distribution (Johnson

1949) to the forestry literature.

The original approach for developing a diameter

distribution model involved (1) estimating the PDF

parameters from sample data using either the method of

moments or maximum likelihood and (2) development of

regression equations to relate the parameter estimates to the

stand variables age, site, and density. These so-called

parameter prediction models (Hyink 1980, Hyink and

Moser 1983) were straightforward to fit. Unfortunately,

functions for relating the PDF parameters to stand

characteristics usually accounted for only a small

percentage of the variation in the parameter estimates (i.e.,

low 2 R values; see Feduccia and others 1979, Smalley and

Bailey 1974). The alternative put forth by Hyink (1980)

was the parameter recovery method. In this method, the

parameters of the distribution function are solved

(recovered) from a system of equations derived by equating

measured or predicted forest stand attributes to their

analytical counterparts from the PDF. Parameter recovery

models are most often based on percentiles of diameter

distribution or moments of diameter distribution. The

parameter recovery method, in general, is superior to the

parameter prediction method for the projection of future

distribution parameters because d.b.h.-frequency

distribution characteristics, such as mean diameter and

diameter variance, can be projected with more confidence

than the distribution parameters themselves. In many

growth-and-yield studies the two methods—parameter

prediction and parameter recovery—are used

simultaneously. That is, some of the PDF parameters are

predicted, and others are solved using a parameter recovery

approach. Over the past two decades there have been many

diameter distribution growth-and-yield systems based on

parameter recovery (possibly in combination with

parameter prediction) using the Weibull distribution (e.g.,

Baldwin and Feduccia 1987, Cao and others 1982, Gove

and Patil 1998, Matney and Sullivan 1982, Zarnoch and

others 1991), and to a lesser extent the SB distribution (e.g.,

Hafley and Buford 1985, Scolforo and Thierschi 1998,

Tham 1988). The goal of this article is to describe the

development and testing of a new alternative approach for

parameter recovery in the SB distribution.

The SB Distribution

Johnson (1949) realized that a random variable X that is

bounded (i.e., has upper and lower limits) could be

transformed to approximate normality by the transformation

[ ] 1 z x x gx γδ ξξλ ln ( ) /( ) ( ) − =+ − +− =

where ln is the natural logarithm and Z   (0,1). Now Z

has PDF

2 1 ( ) exp 2 2

z

n z

π

 −  =    
2

where exp is the base of the natural logarithm. The equation

for Johnson’s SB distribution therefore is

1 1 f x n g x dg x dx ( ) ( ( )) | ( ) / | − − =

or

2

1 exp ln , 2 ( )( ) 2

( ) , 0, , 0, 0

0 otherwise

x

xx x

f x x

δλ ξ γ δ

π ξξ λ ξ λ

ξ ξ λδ γ λ ξ

         −    − +      − +− +−          =  < < + > −∞ < < ∞ > ≥ 







It is characterized by the location parameter ξ, the range

parameter λ, and shape parameters γ and δ.

Hafley and Schreuder (1977) examined the skewness

coefficient, β1 , and kurtosis coefficient, β2, of various

statistical distributions as a measure of the flexibility of the

distributions in regard to their changes in shape. Skewness

is defined as

3/2

1 32 β μμ = /

and kurtosis as

2

2 42 β μμ = /

where

[ ] ( ) ( ) and ( ) r

r μ x E X f x dx f x ∞

−∞ = − ∫ is the PDF of X

Skewness is a departure from symmetry about the mean,

where negative values indicate a distribution with a longer

tail to the left and positive values a longer tail to the right of

the mean. Kurtosis is a measure of the heaviness of the tails

of a distribution; the larger the value of β2 the more

sizeable the tails. There is a myth that kurtosis measures the

peakedness of a density, despite repeated examples to the

contrary (Ali 1974, Johnson and others 1980, Kaplansky

1945). For the standard normal, Z   (0,1), β1 = 0 and

β2 = 3. Moment estimators of β1 and β2 are

3

1 3/2 2

( )

( )

i

i

x x

b n

x x

− =

  −  

∑ ⋅

∑

4

2 2 2

( )

( )

i

i

x x

b

x x

n − =

  −  

∑ ⋅

∑

where n is the sample size.

In figure 1, a graph of the β β 1 2 − space is presented showing

the Weibull, gamma, and lognormal distributions. By

convention the ordinate 2 ( ) β scale is presented upside

down. Certain combinations of β1

 and β2 are

mathematically impossible and occur in the region above

the line β β 2 1 − −=1 0. The three distributions shown are

represented by lines in the β β 1 2 − space. Because the graph

of figure 1 presents β1

, the square of the skewness

coefficient, the positive and negative aspect of a distribution

is not obvious. The lower line of the Weibull plot in figure 1

is generated by negatively skewed shapes. The beta

distribution covers the region between the gamma line, the

impossible region, and the β2 axis. Johnson’s SB distribution

covers the region between the lognormal line, the

impossible region, and the β2 axis. Hence, the beta and SB

distributions cover a broad spectrum of shapes, fitting both

positively and negatively skewed data.

(1)

(2)

β1 - Skewness Squared

01234

β2 - Kurtosis

0

1

2

3

4

5

6

7

8

+

++

+

+

+

+

+

+

+

+

+

+

+

++

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+ +

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

++

+

+

+

+ +

+

+ +

+

+

++

+

+

+

+

+

+

+

+

+

+

+

+

+ +

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+ +

+

++

+

+ + +

+

+

+

+

+

+

+

+ + +

+

+

+

+

+

+

+

+

+

+

+

+

+

+ +

++

+

+

+

+

+

+

+

+

++

+

+

+

+ +

+

+

+

+

+

+

+

+

+ ++ +

+

+

+

+

+

+

+ + +

+

+ +

+

+

++

++ +

+

+

++ +

+

+

+ +

+

+

+

++

+

+

+

+

+ +

+

+

+

+

+

+

+

+

+

+

+

+

+

++

+

+

+

+

+

+

+ +

+

+

+

+

+

+

+

++

+

+

+

+

+

+

+

+

+

+

+

+

+ +

+

++

+

+

+

+

+

+

+

+

+

++

++

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

++

+ +

+

+

+

+

+

+++

+

+ +

+

+

+

+

+

+

+

++

+

+

+

+

++

+

+

+

+

+

+

+ +

+

+

+ +

+

+

+

++

+

+

+

+

+

+

+

+

+ +

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+ +

+

+

+

+

+

+ +

+

+

+

+

+

+

+

+

+

+

+

++

+

+

+

+

+

+ +

+

+

+

+

+

+

+

+

+

+

+

+

+

++

+

+

+

+

+ +

+

+

+

+

+

+

+

+ +

+

+

+

+

+

+

+

+

+

+

++

+

+

+

+

+

+

+ +

+

+

+ +

+ ++

+

+

+ +

+

+

+

+

+

+

+

+

+

+

+

+

+ +

+

+ +

+

+

++

+

+ +

+

+

+

+

+

+

+ +

+

+ +

+

+

+

+

+

+

+

+

+

+

+++

+

+

+

++

+

+

+

+

+

+

+

IMPOSSIBLE REGION

WEIBULL

GAMMA

LOGNORMAL

Figure 1—The β β 1 2 − space. The beta distribution covers the region

between the gamma line and the impossible region. The SB

distribution covers the region between the lognormal line and the

impossible region. The + symbols mark skewness squared-kurtosis

values for the 527 unthinned loblolly pine plot observations.

(3)
3

Because the SB distribution is obtained by a transformation

on a standard normal variate, integration over specific

classes can be accomplished by application of the welltabulated standard normal. Further, the distribution can

easily be extended to multivariate forms (Hafley and Buford

1985, Knoebel and Burkhart 1991). These factors plus the

distribution’s flexibility in the β β 1 2 − space were the salient

points Hafley and Schreuder (1977) used to justify selection

of the SB distribution over other distributions.

Previous SB Models

In previous SB-based models the general approach has been

to specify the distribution minimum and maximum values,

thereby fixing ˆξ = minimum value and ˆ

λ = maximum value

– minimum value, and to recover or otherwise estimate the

two shape parameters. Hafley and Buford (1985) presented

a modeling approach whereby they used a maximum

likelihood type estimating equation for δ, that is,

ˆ ˆ

4 x s

λ δ =

where x s is the standard deviation of x (which is a function

of the first and second moments of a distribution), and for γ

they used

ˆˆ ˆ 2 ˆ ˆ ln ˆˆ ˆ ˆ

m m

m

x x

x

ξλ ξ γ δ λδ λ ξ

−− −   = −  

  + −

where mx is the mode of x. Other researchers have relied on

the two-percentile method to recover γ and δ. In Newberry

and Burk (1985) and Knoebel and Burkhart (1991) the 50th

and 95th percentiles were used. From the formulation of the

SB distribution,

95 50

95 50

95 50

ln and ln X X Z Z

X X

ξ ξ γδ γδ ξλ ξλ

  − − =+ =+     +− +−

where Z95 and Z50 represent the standard normal values

corresponding to the cumulative frequencies of 95 percent

and 50 percent of the standard normal distribution, and X95

and X50 are the 95th and 50th percentiles of the observed

distribution. Because Z50 = 0, we have

50

50

ˆ ˆ ˆ ln ˆ ˆ

X

X

ξ γ δ ξ λ

  − = −  

  + −

Substituting into the equation for Z95 and simplifying yields

95 50

95 50

1.645 ˆ ˆ ˆ ln ln ˆˆ ˆˆ

X X

X X

δ

ξ ξ

ξλ ξλ

=    − −    −

   +− +−

Newbury and Burk (1985) correctly pointed out that these

two percentile points, the 50th and 95th, do not necessarily

provide the optimal estimates for γ and δ.

Of course, many approaches have been used to determine

the SB parameters, such as the four-percentile points

method, linear and nonlinear regression methods, and

maximum likelihood. These have been reviewed and

compared by Zhou and McTague (1996) and Kamziah and

others. (1999). The purpose of my brief review of the above

methodology is to provide a contrast against the

methodology presented in the next section.

Development of the SB Parameter

Recovery System

The parameter recovery method ensures compatibility

between stand characteristics generated from a distribution

function and predicted from regressions. To apply the

parameter recovery method, a system of equations involving

certain tree stand characteristics and the SB parameters must

be developed. If we define a new variable Y with variate

values

y x = − ( )/ ξ λ

then the SB distribution [equation (1)] can be expressed as

2

1 ( ) exp 0.5 ln , 0 1

2 (1 ) 1

y fy y

yy y

δ γ δ

π

        = − + <<       − −        

where

ln[ /(1 )] y z yy =+ − γ δ

Solving equation (5) for y gives

( )/ 1 (1 ) ( ) z y e gz −− − γ δ =+ =

Setting z = 0 in equation (6) results in the median value of

Y: / 1

median y e (1 ) γ δ − = +

hence

median γ δ = − ln(1/ 1) y

Because Y is a function of Z [see equation (6)],

( ) [ ( )] r r EY E gZ = . As Z is a unit normal variable, the rth

noncentral moment of Y is

1 2 ( ) exp( / 2) /[1 exp(( ) / )]

2

r

r μ γδ y z z dz

π

∞

−∞

′ = − +− ∫

(4)

(5)

(6)

(7)

(8)
4

which can be evaluated using a formula due to Goodwin

(1949). Hill and others (1976) give an excellent algorithm

for evaluating SB moments based on Draper’s (1952) form

of Goodwin’s integral.

Let us proceed by solving equation (4) for x, which gives

x y = + λ ξ . Taking the statistical expectation [i.e.,

EX E Y () ( ) = + λ ξ ] we obtain

1 x y = + λμ ξ ′( )

Quadratic mean d.b.h. q ( ) D , number of trees per unit area

(N), and basal area per unit area (B) are related by the

equation (see Clutter and others 1983, p. 89)

q D BN = /( ) , thus 2 B ND = ⋅ q where is for units

conversion from d.b.h.2

 to basal area (= π / 40,000 for

converting d.b.h.2

 in cm2 to area in m2) and 2 2

q D EX = ( ). In

terms of our new variable Y, 2 2 EX E Y ( ) [( ) ] = += λ ξ 2 2

2 1 λ μ ξλμ ξ ′ ′ () 2 () y y + + . Now B, expressed as a function

of the first and second moments of Y, becomes

2 2

2 1 BN y y = ++  [ () 2 () ] λ μ ξλμ ξ ′ ′

The relationship in equation (7) is used to eliminate γ in

equations (9) and (10) by substitution into equation (8). We

are then left with a system of two nonlinear equations and

two unknowns. Given estimates of B, N, median tree

diameter median ( ) d , average tree diameter mean ( ) d , and the ξ

parameter, plus the relationships in equations (9) and (10),

we can iteratively solve for parameters λ and δ. The

parameter γ is then determined from equation (7).

An attempt was made to define three equations in three

unknowns so that the ξ parameter could be recovered as

well. Taking the variance of equation (4) gives

2 2 var( ) ( ) ( ) X yy λμ μ 2 1 = −   ′ ′  

However, this function is a combination of equations (9)

and (10), which is readily seen by subtracting the square of

equation (9) from equation (10). A system based on these

three equations cannot be solved. Next kurtosis was tried,

which, expressed in terms of noncentral moments, is

2 4

4 13 1 2 1

2 2 2

2 1

() 4 () () 6 () () 3 ()

() ()

y yy y y y

y y

μ μμ μ μ μ β

μ μ

′ ′′ ′ ′ ′ −+ − =   ′ ′ −  

However, the resulting complexity of this system of three

nonlinear equations posed such convergence problems that

it was abandoned.

Undoubtedly overall stand attributes such as volume and

basal area are less sensitive to errors in ξ than to errors in

the other three parameters. Hence, specifying ξ and

recovering the remaining parameters is probably the most

expedient approach.

Solution Technique

The solution of a system of nonlinear equations is a difficult

problem. A two-dimensional Newton-Raphson procedure is

relatively straightforward to program and has quadratic

convergence. The two-dimensional Newton-Raphson

method can be written in vector notation as (Ralston and

Rabinowitz 1978, p. 360)

2 1

(2) (2) 1

1

2 1 2

(1) (1)

1

i

i i

f f

f x x

g f f f

x x

+

=

∂ ∂    −     ∂ ∂   = −    ∂ ∂    −       ∂ ∂ x x

x x

where the determinant g is

1 1

(1) (2)

2 2

(1) (2)

f f

x x

g f f

x x

∂ ∂

∂ ∂ = ∂ ∂

∂ ∂

The key to using equation (11) is correctly specifying the

partial derivatives of 1f and 2f with respect to (1) x (the λ

parameter) and (2) x (the δ parameter). I shall denote the

relationship in equation (9) as 1f and in equation (10) as 2f .

The partial derivatives of equations (9) and (10) depend on

the partial derivatives of the moments of Y. Let us define the

pseudo-moment ( ) r μ′ y as

2 1/ 2 exp( / 2) /[1 exp(( ) / )]r π γδ z z z dz

∞

−∞

− +− ∫

The partial derivatives of the moments of Y (see appendix)

can be shown to be functions of ( ) r μ′ y and ( ) r μ′ y . Using the

results of equations (A1) and (A2), the partial derivatives of

1f and 2f are:

1 ( ) 1 2

1

median

() () ( ) f y y

y d

λμ μ μ λ λξ

∂ ′ ′ − = − ′ ∂ −+

1 ( ) 1 2

2

f λμ μ () () y y

δ δ

∂ − − ′ ′ = ∂

 

(9)

(10)

(11)

(12)

(13)

(14)
5

( )( )

( )

23 12

2

median

2 1

2 () () () ()

 2 ( ) ( )

yy yy

f d N

y y

λλμ μ ξμ μ

λ ξ λ

λμ ξμ

  − −+ −   ′′ ′′   ∂   =   − +

∂   + + ′ ′  

⋅

( )( ) 23 12 2

2

f 2 () () () () yy yy N λλμ μ ξμ μ

δ δ

∂   − −+ −   ′′ ′′     =   ∂    

⋅    

Substituting equations (13) through (16) into equation (11)

gives us the two-dimensional Newton-Raphson iteration

formula to solve the system of equations. The L1 norm

2

1

1

|| || i

i

f =

f = ∑

was used to define convergence as

1 || || 0.002 f <

Data

Data came from 287 plots (varying in size from 0.04 to

0.3 ha) established in unthinned loblolly pine (Pinus taeda

L.) plantations in north and central Louisiana, southern

Arkansas, southeast Texas, and southwest Mississippi.

Periodic remeasurements of stand conditions were taken

every 4 to 6 years on 125 of the plots. The database

consisted of 527 stand observations, 240 of which were

growth observations (remeasurements). Diameter at breast

height (1.37 m) to the nearest 0.25 cm was measured for

each tree on the plot. The average height of dominant and

codominant trees was obtained on each plot at each

measurement to determine site index. Stand ages ranged

from 5 to 45 years, trees per hectare ranged from 119 to

3,039, basal area per hectare varied from 1.15 to 48.2 m2

,

and average height of dominants and codominants varied

from 5.2 to 27.4 m.

Figure 1 presents a plot of the 527 ( 1 b , 2 b ) values, calculated

from equations (2) and (3) on the diameter measurement

values. Few points fall near the line defined by the Weibull

distribution. Clearly, the SB distribution covers most of the

( 1 b , 2 b ) values and therefore should be appropriate for

modeling these stand diameter distributions.

The 0th Percentile for the Distribution

Minimum Value

Estimates of the location parameters of the recovered

distributions are needed. In previous studies (Cao and

others 1982, Matney and Sullivan 1982, Zöhrer 1972)

regressions were developed on the sample plot minimum

diameters, then subjective adjustments, such as dividing the

prediction by 2, were made to estimate the threshold

parameters. Other studies (Dell and others 1979, Feduccia

and others 1979, Smalley and Bailey 1974) developed

direct estimation of threshold parameters from regressions

on maximum likelihood estimates (MLEs). Objective

methods are preferable, so maximum likelihood was used to

solve for the 527 ξ parameters. However, the MLEs were

unsatisfactory. Values of ˆξ were nearly 0 in most cases and

often were in conflict with the natural progression of

diameter over time.

Looking at plots of the empirical CDFs, I observed that the

curves tended to straighten at the ends. I took the three

smallest plot diameter values and linearly regressed down to

the 0th percentile value (fig. 2). Percentiles were calculated

on the ordered d.b.h. values d1, d2, . . ., dn as

( ) /( 1) 100 pd i n i = +× . This technique gave reasonable

diameter values for the minimum distribution values.

(15)

(16)

DBH (cm)

8 10 12 14 16 18 20 22 24 26 28 30

Percentile

0

10

20

30

40

50

60

70

80

90

100

x x x x x

x

x

x

x

x

x

x

x

x

x

x

x

x

x x

x x

x

x

Figure 2—Empirical CDF of a loblolly pine study plot with a linear

regression line fitted to the three smallest values showing extrapolation

to the 0th percentile value.
6

Stand Equations

Compatible projection models for number of trees surviving

and basal area were fitted to the sample data. Additionally,

prediction equations for median d.b.h., mean d.b.h., and the

ξ parameter also were fitted. Table 1 shows the equations

that predict stand attributes. The loblolly pine site index

curve used was that of Popham and others (1979).

In addition to stand attribute values, starting values are

required to begin the iteration of equation (11). For the

range parameter, the following equation gave good starting

values:

2 2 initial 12.29 0.2562 0.0123 ( / ) 0.00338 / x x λ = + ⋅− ⋅ + ⋅⋅ s NA s NA

where 2

x s is the variance of d.b.h. in cm2

, N is number of

trees per hectare, and A is stand age in years. For δ, the

value 1.2 was determined to be a good starting point, based

on the observation that δ often occurs between 1 and 1.5 for

SB diameter distributions.

Results

A FORTRAN program was written to facilitate use of the

parameter recovery model and to provide a means of

Table 1—Equations for predicting unthinned loblolly plantation stand attributes

Equationa

6 1.06412 2.38353 2.47812 2

2 1 2 1 11 1 1

b

ˆ 1.05409 10 ( ) ( / ) exp(4.56524 10 )

240; 0.964; RMSE 83.5

N N A A N B HD B

n FI

− − =− −

== =

⋅ ⋅ ⋅ ⋅ ⋅⋅

2 0.42213 4.38106 10 2.49276 21 12.69457

2 1 211 1 2 1

ˆ 405.39237( ) exp( 5.10433 10 )

240; 0.942; RMSE 2.3674

B B A A B HD A B

n FI

− ⋅ − − =+ − −

== =

⋅ ⋅ ⋅ ⋅ ⋅⋅

2

median q

ˆ 0.51862 1.02175 1.60392 10 ( 1 if 20 else 0)

527; 0.988; RMSE 0.5428

d D AI I A I

n FI

− =− + − = ≥ =

== =

⋅ ⋅⋅

q mean ln( ) 1.95405 0.43658 ln( )

527; 0.999; RMSE 0.1732

D d HD

n FI

− =− +

== =

2 22

q q mean ˆ 0.90155 0.018167 0.27733( ) 0.43026

523; 0.757; RMSE 2.1996

D D d HD

n FI

ξ =− + − − +

== =

a Notation: A = stand age in years, B = basal area in m2

 per hectare, dmean = average d.b.h. in cm, dmedian =

median d.b.h. in cm, Dq = quadratic mean d.b.h. in cm, HD = average height of dominants and codominants

in m, N = number of surviving trees per hectare, and ξ = distribution minimum d.b.h. in cm. Subscripts refer

to time of measurement.

b Statistics: fit index, 2 2 1 ( )/ ( ) ˆ FI =− − − ∑ ∑ yy yy ii ii ; root mean square error, RMSE = 2 ( ) /( ) ˆ i i ∑ y y np − − ;

where n = number of observations, p = number of parameters, yi

 = observed value, y = observed mean value,

and ˆi y = predicted value.

b
7

checking the model against the data. Inputs into the

program were initial age, initial basal area, initial number of

trees, initial height of dominant and codominant trees or

base age 25 site index, the number of projection periods,

and the time length between projection periods. The

program provided the recovered parameter values, current

and future stand tables, and histograms.

Initial stand conditions for each of the 287 plots were input

into the program, and 287 current and 240 future stand

tables were generated and compared against the observed

distributions. A convergent solution was obtained on all 527

observations. On the average, convergence was obtained in

3 to 4 iterations. From a purely mechanical viewpoint, the

SB parameter recovery model worked extremely well. The

degree of correspondence between the actual and predicted

d.b.h.-frequency distributions was evaluated next.

For all observations, differences were calculated between

observed and predicted number of trees by 2.5-cm d.b.h.

classes. The generated SB distributions accounted for 77

percent of the observed variation in number of trees by

d.b.h. class. The standard deviation was 67 trees per

hectare. Residual analysis showed no apparent prediction

bias except in the smallest diameter class, where numbers

were slightly underestimated.

The 2 χ goodness-of-fit test, using 2.5-cm classes, was

performed on each of the observations. These tests resulted

in rejection of the hypothesis that diameter distribution is

from an SB distribution in 56 cases at the α = 0.01 level.

While this number is considerably larger than the expected

number of rejections (0.01 527 5) × ≅ , it was not

unanticipated. The graph of skewness squared-kurtosis

values in figure 1 clearly shows some values outside the

range covered by the SB distribution. Of the 56 cases

rejected, 50 had values [i.e., ( 1 b , 2 b ) points] that fell below

the lognormal line. Taking these 50 observed plot

distributions into account (56 50 6) − = , the 2 χ test results

are very reasonable. It was concluded, therefore, the

recovered SB distributions provided good approximations of

the observed diameter distributions.

Discussion

A parameter recovery framework was developed using the

stand level values of mean d.b.h. and basal area, in

conjunction with median d.b.h., and the first and second

moments of the SB distribution. The resulting system of two

nonlinear equations [equations (9) and (10)] plus equation

(7), which form the parameter recovery model, allowed for

the determination (i.e., recovery) of the range and both

shape parameters of the SB distribution. The approach of

Hafley and Buford (1985) and the two-percentile method

used by Newbury and Burk (1985) and Knoebel and

Burkhart (1991) do not require iterative procedures.

However, Kamziah and others (1999) were critical of the

two-percentile method because “. . . the parameters γ and δ

were estimated based only on two percentiles, and λ was

not estimated simultaneously from these percentiles.” They

felt it was important to interconnect the value of the range

parameter to the overall shape of the distribution, because

the range and shape parameters provide feedback to each

other (a gestalt, if you will) and, therefore, should be

estimated simultaneously.

This same criticism can be applied to the approach of

Hafley and Buford (1985). This new alternate methodology

provides that link between the range and shape parameters.

The methodology can be extended to the SBB distribution, a

bivariate distribution that is commonly used for describing

simultaneously diameter and height structure (Tewari and

Gadow 1999). In that case, the marginal distributions for

diameter and height would be determined as outlined, but a

strategy for determining the correlation parameter to tie the

marginal distributions together would have to be developed.

The wealth of literature on diameter distribution modeling

suggests no one system will universally fit all situations. In

this study approximately 10 percent of the plots could not

adequately be modeled using the SB distribution. A better

understanding of the biological relationships that

underscore stand development over time is necessary to

further develop existing approaches. Nevertheless, the SB

parameter recovery system developed in this paper is

tractable for use in stand modeling and yielded good

representations of unthinned loblolly pine stands.

Literature Cited

Ali, M.M. 1974. Stochastic ordering and kurtosis measure. Journal of the

American Statistical Association. 69: 543-545.

Bailey, R.L.; Dell, T.R. 1973. Quantifying diameter distributions with the

Weibull function. Forest Science. 19(2): 97-104.

Baldwin, Jr., V.C.; Feduccia, D.P. 1987. Loblolly pine growth and yield

prediction for managed west gulf plantations. Res. Pap. SO-236. New

Orleans: U.S. Department of Agriculture, Forest Service, Southern

Forest Experiment Station. 27 p.

Cao, Q.V.; Burkhart, H.E.; Lemin, Jr., R.C. 1982. Diameter

distributions and yields of thinned loblolly pine plantations. FWS-1-82.

Blacksburg, VA: Virginia Polytechnic Institute and State University.

School of Forestry and Wildlife Resources. 62 p.
8

Clutter, J.L.; Bennett, F.A. 1965. Diameter distribution in old-field slash

pine plantations. Rep. 13. Macon, GA: Georgia Forest Research

Council. 9 p.

Clutter, J.L.; Fortson, J.C.; Pienaar, L.V. [and others]. 1983. Timber

management: a quantitative approach. New York: John Wiley. 333 p.

Dell, T.R.; Feduccia, D.P.; Campbell, T.E. [and others]. 1979. Yield of

unthinned slash pine plantations on cutover sites in the West Gulf

region. Res. Pap. SO-147. New Orleans: U.S. Department of

Agriculture, Forest Service, Southern Forest Experiment Station. 84 p.

Draper, J. 1952. Properties of distributions resulting from certain simple

transformations of the normal distribution. Biometrika. 39: 290-301.

Feduccia, D.P.; Dell, T.R.; Mann, Jr., W.F. [and others]. 1979. Yields of

unthinned loblolly pine plantations on cutover sites in the west gulf

region. Res. Pap. SO-148. U.S. Department of Agriculture, Forest

Service, Southern Forest Experiment Station. 88 p.

Goodwin, E.T. 1949. The evaluation of integrals of the form

2

( ) x f x e dx ∞ −

∫−∞ . Proceedings of the Cambridge Philosophical Society.

45: 241-245.

Gove, J.H.; Patil, G.P. 1998. Modeling the basal area-size distribution of

forest stands: A compatible approach. Forest Science. 44(2): 285-297.

Hafley, W.L.; Buford, M.A. 1985. A bivariate model for growth and yield

prediction. Forest Science. 31(1): 237-247.

Hafley, W.L.; Schreuder, H.T. 1977. Statistical distributions for fitting

diameter and height data in even-aged stands. Canadian Journal of

Forest Research. 7: 481-487.

Hill, I.D.; Hill, R.; Holder, R.L. 1976. Algorithm AS 99. Fitting Johnson

curves by moments. Applied Statistics. 25: 180-189.

Hyink, D.M. 1980. Diameter distribution approaches to growth and yield

modeling. In: Brown, K.M.; Clarke, F.R., eds. Forecasting forest stand

dynamics. Thunderbay, Ontario: Lakehead University, School of

Forestry: 138-163.

Hyink, D.M.; Moser, Jr., J.W. 1983. A generalized framework for

projecting forest yield and stand structure using diameter distributions.

Forest Science. 29(1): 85-95.

Johnson, M.E.; Tietjen, G.L.; Beckman, R.J. 1980. A new family of

probability distributions with applications to Monte Carlo studies.

Journal of the American Statistical Association. 75: 276-279.

Johnson, N.L. 1949. Systems of frequency curves generated by methods

of translation. Biometrika. 36: 149-176.

Kamziah, A.K.; Ahmad, M.I.; Lapongan, J. 1999. Nonlinear regression

approach to estimating Johnson SB parameters for diameter data.

Canadian Journal of Forestry. 29: 310-314.

Kaplansky, I. 1945. A common error concerning kurtosis. Journal of the

American Statistical Association. 40: 259.

Knoebel, B.R.; Burkhart, H.E. 1991. A bivariate distribution approach to

modeling forest diameter distributions at two points in time.

Biometrics. 47(1): 241-253.

Matney, T.G.; Sullivan, A.D. 1982. Compatible stand and stock tables for

thinned and unthinned loblolly pine stands. Forest Science. 28(1):

161-171.

Newbury, J.D.; Burk, T.E. 1985. SB distribution-based models for

individual tree merchantable volume-total volume ratio. Forest Science.

31: 389-398.

Popham, T.W.; Feduccia, D.P.; Dell, T.R. [and others]. 1979. Site index

for loblolly plantations on cutover sites in the west gulf coastal plain.

Res. Note SO-250. New Orleans: U.S. Department of Agriculture,

Forest Service, Southern Forest Experiment Station. 7 p.

Ralston, A.; Rabinowitz, P. 1978. A first course in numerical analysis.

2nd ed. New York: McGraw-Hill. 556 p.

Scolforo, J.R.S.; Thierschi, A. 1998. Estimativas e testes da distribuição

de frequência diâmétrica para Eucalyptus camaldulensis, através da

distribuição SB, por diferentes métodos de ajuste. Scientia Forestalis.

54: 93-106.

Smalley, G.W.; Bailey, R.L. 1974. Yield tables and stand structure for

loblolly pine plantations in Tennessee, Alabama, and Georgia

highlands. Res. Pap. SO-96. New Orleans: U.S. Department of

Agriculture, Forest Service, Southern Forest Experiment Station. 81 p.

Tewari, V.P.; Gadow, K.V. 1999. Modelling the relationship between tree

diameters and heights using SBB distribution. Forest Ecology and

Management. 119: 171-176.

Tham, Å. 1988. Estimate and test frequency distributions with the

Johnson SB function from stand parameters in young mixed stands after

different thinning treatments. In: Ek, A.R.; Shifley, S.R.; Burk, T.E.,

eds. Forest growth modelling and prediction: Proceedings of the

IUFRO Conference. Gen. Tech. Rep. NC-120. Minneapolis: U.S.

Department of Agriculture, Forest Service, North Central Forest

Experiment Station: 255-262. Vol. 1.

Zarnoch, S.J.; Feduccia, D.P.; Baldwin, Jr., V.C.; Dell, T.R. 1991.

Growth and yield predictions for thinned and unthinned slash pine

plantations on cutover sites in the west gulf region. Res. Pap. SO-264.

New Orleans: U.S. Department of Agriculture, Forest Service, Southern

Forest Experiment Station. 32 p.

Zhou, B.; McTague, J.P. 1996. Comparison and evaluation of five

methods of estimation of the Johnson system of parameters. Canadian

Journal of Forestry. 26: 928-935.

Zöhrer, F. 1972. The beta-distribution for best fit of stem-diameterdistribution. In: 3rd Conference advisory group of Forest Statisticians:

Proceeding of the IUFRO Conference. INRA Publ. 72-3. Paris: Institut

National Recherche Agronomique: 91-106.
9

Appendix

The partial derivatives of the moments of Y [ ( ) r μ′ y , see equation (8)] with respect to λ and δ are derived below. These

derivatives are needed for the computation of the partial derivatives of equations (9) and(10). First, performing a change-ofvariable from median y to dmedian [where we let d, for tree diameter, replace x in equation (4)], we have from equation (7) that

[ ] median γδ λ ξ = −− ln /( ) 1 d . Substituting for γ where necessary, we have

[ ]

2

median

1

2

1

median

1

median

( ) exp( / 2) ln[ /( ) 1]

2 [1 exp(( ) / )]

1 exp( / 2)

2 [1 exp(( ) / )]

() ()

exp(( ) / )

exp(( ) / )

r

r

r

r r

yr z d z D dz

z

r z dz

d z

r y y d

z

z

λ

μ δλ ξ

λ γδ δ π

λ ξ γδ π

μ μ λ ξ

γ δ

γ δ

∞

+

−∞

∞

+

−∞

+

∂ − ′   − −− = −   ∂ +−  

− = − −+ + −

=− − ′ ′ − +

⋅ − ⋅

⋅ −

⋅

∫

∫

and

2

median

1

2

2 1

( ) exp( / 2) ln[ /( ) 1]

2 [1 exp(( ) / )]

1 exp( / 2)

2 [1 exp(( ) / )]

exp(( ) / )

exp(( ) / )

r

r

r

yr z d z D dz

z

r z z dz

z

z

z

δ

μ δλ ξ

δ γδ δ π

δ γδ π

γ δ

γ δ

∞

+

−∞

∞

+

−∞

∂ − ′   − −− = −   ∂ +−  

− = − + −

⋅ − ⋅

⋅ − ⋅

∫

∫

Because of the extra z, the integral does not reduce to the difference of two SB moments as it did in the derivative with respect

to λ. Using the pseudo-moment ( ) r μ′ y defined in equation (12), we can express the partial derivative as

2 [ ] 1

( ) () () r

r r

y r y y μ μ μ δ δ +

∂ ′ =− − ′ ′ ∂ ⋅  

(A1)

(A2)
10
Parresol, Bernard R. 2003. Recovering parameters of Johnson’s SB distribution. Res. Pap. SRS-31.

Asheville, NC: U.S. Department of Agriculture, Forest Service, Southern Research Station. 9 p.

A new parameter recovery model for Johnson’s SB distribution is developed. This latest alternative approach

permits recovery of the range and both shape parameters. Previous models recovered only the two shape

parameters. Also, a simple procedure for estimating the distribution minimum from sample values is presented.

The new methodology employs the median and the first and second moments of the distribution. The methodology is demonstrated by modeling diameter distributions of unthinned loblolly pine plantations. Compatible

equations for projecting per-hectare values of number of trees and basal area from initial stand conditions are

presented, as well as equations for predicting median diameter, mean diameter, and the location parameter.

Given estimates of these five stand attributes, the range and the two shape parameters of the SB distribution can

be recovered. The 2 χ goodness-of-fit test rejected 56 cases out of 527 for conformance to an SB distribution.

Though the SB distribution is very flexible in terms of distribution shape, about 10 percent of the loblolly

plantation observations did not follow this distribution. Nonetheless, deviation analysis showed reasonable

results, with 77 percent of the variation explained in current and projected distributions (numbers of trees by

2.5-cm diameter class). Overall, the recovered SB distributions provided good approximations of the observed

diameter distributions.

Keywords: Diameter distributions, kurtosis, Newton-Raphson method, moments, parameter recovery model,

percentiles, Pinus taeda, skewness.
The Forest Service, United States Department of

Agriculture (USDA), is dedicated to the principle of

multiple use management of the Nation’s forest resources

for sustained yields of wood, water, forage, wildlife, and recreation.

Through forestry research, cooperation with the States and private

forest owners, and management of the National Forests and National

Grasslands, it strives—as directed by Congress—to provide

increasingly greater service to a growing Nation.

The USDA prohibits discrimination in all its programs and activities

on the basis of race, color, national origin, sex, religion, age,

disability, political beliefs, sexual orientation, or marital or family

status. (Not all prohibited bases apply to all programs.) Persons

with disabilities who require alternative means for communication of

program information (Braille, large print, audiotape, etc.) should

contact USDA's TARGET Center at (202) 720-2600 (voice and

TDD).

To file a complaint of discrimination, write USDA, Director, Office

of Civil Rights, Room 326-W, Whitten Building, 1400

Independence Avenue, SW, Washington, D.C. 20250-9410 or call

(202) 720-5964 (voice and TDD). USDA is an equal opportunity

provider and employer.