A comparison of procedures for estimating the parent probability distribution from a given set of fractiles

European Journal of Operational Research 120 (2000) 657±670

www.elsevier.com/locate/orms

Theory and Methodology

A comparison of procedures for estimating the parent probability distribution from a given set of fractiles Hon-Shiang Lau a

a,*

, Amy Hing-Ling Lau a, John F. Kottas

b,1

College of Business Administration, Oklahoma State University, Stillwater, OK 74078, USA b School of Business, College of William and Mary, Williamsburg, VA 23185, USA Received 1 February 1998; accepted 1 September 1998

Abstract On one hand, eliciting subjective probabilities (fractiles) is a well-established procedure. On the other hand, knowledge of a subjective variable's central moments or distribution function is often assumed. However, the problem of converting elicited fractiles into the required moments or distribution function has been largely ignored. We show that this conversion problem is far from trivial, and that the most commonly used conversion procedures often produce huge errors. Alternative procedures are proposed; the ``Tocher's curve'' and ``linear function of fractiles'' methods are shown to be much more accurate than the commonly used procedures. Ó 2000 Elsevier Science B.V. All rights reserved. Keywords: Risk analysis; Statistical distribution theory; Forecasting

1. Introduction On one hand, a huge literature exists on how one should elicit the fractiles of a subjective probability distribution (see, e.g., reference lists in van Lenthe, 1994; Lau et al., 1996). On the other hand, since Hertz's (Hertz, 1964) classical paper, many risk/decision models now assume the knowledge of the moments and/or the distribution functions of their subjectively estimated variables (e.g., the future annual cash ¯ows of a new project). However, little has been said about how one should convert the elicited fractiles into the moments and/or distribution functions required by the decision models. The implicit assumption is that the conversion problem is trivial. For example, Hertz and Thomas (1983, p. 161) stated that the cumulative distribution function can be obtained simply as ``a smooth curve drawn by eye through the plotted (fractile) points''. This assumption is further reinforced in recent years by the increasingly popular ``DecisionTools

* 1

Corresponding author. Tel.: 405-744-5105; fax: 405-744-5180; e-mail: [email protected] Authors' names listed in reverse alphabetical order.

0377-2217/00/$ - see front matter Ó 2000 Elsevier Science B.V. All rights reserved. PII: S 0 3 7 7 - 2 2 1 7 ( 9 9 ) 0 0 0 0 2 - 8

658

H.-S. Lau et al. / European Journal of Operational Research 120 (2000) 657±670

Table 1 Average moment deviations from ®ve ®tting methods (2000 random beta test distributions) Fitting method Piecewise linear approximation

BESTFIT-®tted distribution

For {S4 } ˆ {2, 3, 4, 5, 6, 7, 8} Dl (% of r) Dr (% of r) Db1 Db2

162 793 13.7 14.7

162 793 1.16 9.14

For {S3 } ˆ {1, 3, 4, 5, 6, 7, 9} Dl (% of r) Dr (% of r) Db1 Db2

36 342 42.1 51.6

36 342

CPU time (seconds) a b

0.0002

a a

b b

0.0002 + ?

Akima's spline

Tocher's curve

Linear function of fractiles

1.4 4.9 0.67 1.67

0.18 0.82 0.120 0.310

0.021 0.39 0.126 0.143

0.9 2.7 0.10 0.47

0.04 0.13 0.048 0.163

0.017 0.36 0.123 0.103

0.002

0.02

.0

20 test problems for BESTFIT-®tted b1f and b2f . Did not investigate.

Suite'' software package (Palisade Corporation, 1996). The BESTFIT component of this ``suite'' produces for any given set of fractiles a ®tted distribution function, which then serves as an input to their risk simulation package. The purposes of this paper are to: (i) show that the conversion problem is far from trivial and ``solved'', and well-known procedures such as the BESTFIT software often produce huge errors; (ii) brie¯y describe various alternative procedures of conversion; (iii) show that these alternatives di€er considerably in accuracies; (iv) identify the most promising alternatives for adoption and further improvement. Section 2 will explain our approach for evaluating any given conversion procedure. Only procedures that can be rapidly (say, < 10 seconds) implemented by a PC are considered in this study; Sections 3±6 describe brie¯y each of these conversion procedures. Section 7 reports the procedures' performance; Table 1 below serves as a convenient preview of the considerable di€erences among the procedures' accuracies. Section 8 discusses future investigation directions. 1.1. De®nition of symbols k

For a random variable x, de®ne l ˆ E…x† ˆ x's expected value; lk ˆ E…x ÿ l† ˆ x's kth central moment p …k > 1†; r ˆ l2 ; a1 ˆ l3 =r3 ˆ x's skewness measure; b1 ˆ a21 ; b2 ˆ l4 =r4 ˆ x's kurtosis measure. For our purpose, the term ``®rst four moments'' (``€m'') refers to collectively x's l, r, b1 and b2 . De®ne also f …x†, F …x† and F ÿ1 …p† as, respectively, x's density function, cumulative distribution function (``cdf'') and inverse cdf. Let xa denote x's ath fractile; thus, x0 and x1 denote respectively x's minimum and maximum points. 2. Approach for evaluating fractile-to-distribution conversion procedures A large number of known ``test distributions'' are selected. For each test distribution, its fractiles are determined. These fractiles become inputs to the conversion procedure to be tested, and one then evaluates how close the conversion-procedure-®tted distribution is to the known test-distribution. Three basic components of this approach are discussed below.

H.-S. Lau et al. / European Journal of Operational Research 120 (2000) 657±670

659

2.1. Selection of test distributions A conversion-procedure's performance should be tested with distributions of various shapes. A standard measure of the ``shape variety'' of a group of distributions is their (b1 , b2 ) coverage (see, e.g., Hahn and Shapiro, 1967). Distributions such as the normal, exponential, gamma and Weibull cover only a point or a line of (b1 , b2 )s (see Fig. 1). The four-parameter Pearson Type-I (commonly known as ``beta'') distribution f …x† ˆ

x ÿ U pÿ1
…V ÿ x†qÿ1 B…p; q†
…V ÿ U †

p‡q‡1

;

U PxPV ;

p 6 0; q 6 0

…1†

(where B(.,.) is the beta function) covers a (b1 , b2 )-area (instead of a point or line), and is widely perceived in the management science literature as one that can represent a ``suciently wide'' variety of shapes. J- or Ushaped beta distributions have either p 6 1 and/or q 6 1 in Eq. (1), and ``bell-shaped'' beta distributions have p > 1 and q > 1 (see Johnson et al., 1995). Note that ``bell-shaped'' is used loosely here to mean an unimodal ``inverted-U'' that is often asymmetrical. It is generally recognized that J- and U-shaped distributions are much less relevant than bell-shaped ones in modelling subjective distributions. Therefore this study reports only results for bell-shaped beta distributions; in other words, we will assume ``a world of bell-shaped beta distributed random variables''. We have, however, replicated all experiments using the entire range of beta f …x†s and veri®ed that they lead to the same conclusions. To simplify the evaluation of the conversion procedure's performance (clari®ed in the subsection after next), we will use standardized beta distributions with l ˆ r ˆ 1. Each conversion procedure will be tested with 2000 beta distributions, each distribution corresponding to a di€erent (b1 , b2 ) that is randomly

Fig. 1. (b1 , b2 ) diagram (adapted from Hahn and Shapiro, 1967).

660

H.-S. Lau et al. / European Journal of Operational Research 120 (2000) 657±670

generated from the (b1 , b2 )-area of the bell-shape beta (see Fig. 1). For example, with the ®rst randomly generated (b1 , b2 ) ˆ (1.56, 5.095), x's (the test distribution's) ``true'' €m are: l ˆ 1; r ˆ 1; b1 ˆ 1:560; b2 ˆ 5:095:

…2†

Substituting these €m into standard closed-form formulas (e.g., from Johnson et al., 1995) gives the betadistribution's parameters (see Eq. (1)) as …U ; V ; p; q† ˆ …ÿ0:451; 23:085; 1:913; 29:127†:

…3†

Note that we need only to use distributions with di€erent (b1 , b2 )s but not with di€erent (l)s and (r)s because it can be easily shown that a conversion procedure's accuracy is a€ected by a distribution's shape but not by its location or scale. Similarly, a conversion-procedure's accuracy remains unchanged when a positively skewed distribution is to its negatively skewed mirror image. Hence we will only conpswitched  sider distributions with a1 ˆ ‡ b1 . 2.2. Selection of fractile-levels to be considered A conversion procedure's performance may be a€ected by the particular fractile-levels (i.e., the ``a''s) used in the fractile-set {xa s} provided. Therefore, the conversion-procedures' performance should be investigated with respect to a suciently wide variety of ``likely and adequate'' fractile sets. An increasingly large number of theoretical, empirical and case studies (e.g., Selvidge, 1980; Alpert and Rai€a, 1982; Solomon, 1982; Lau et al., 1996) support the notion that a subjective distribution should be constructed by eliciting its median (i.e., x0:5 ) plus 4 or 6 symmetrical fractiles from the following set of 9 fractile-levels (ai )s: Index from fSg ˆ fi; i ˆ 1; . . . ; 9g Fractiles Level ai

1 0:01

2 0:05

3 0:10

4 0:25

5 0:50

6 0:75

7 0:90

8 9 0:95 0:99

…4†

From the index set fSg ˆ fi; i ˆ 1; . . . ; 9g, one can form 10 di€erent subsets with 4 or 6 symmetrical (ai )s (plus ai ˆ 0.5), they are designated ffSk g; k ˆ 1 to 10g and enumerated in Appendix A. All 10 fSk gs are considered in our study (there is no consensus in the literature on which fSk g is the best), but for brevity only fS4 g de®ned in (5a) and (5b) will be used in Sections 2±6 as illustration. For fS4 g, the corresponding fractiles given in (5c) for the beta f …x† de®ned by Eq. (2) or Eq. (3) can be computed with the IMSL (1987) subroutine BETIN: Set {S4 }, with i ˆ Fractile level ai Fractile value x…i† ˆ xai

2 0.05 ÿ0.198

3 0.10 ÿ0.066

4 0.25 0.252

5 0.50 0.783

6 0.75 1.515

7 0.90 2.355

8 0.95 2.941

(5a) (5b) (5c)

2.3. Evaluating the closeness between a test distribution and its ®tted distribution The fractiles in (5) for a given test distribution ft …x† are the inputs to the conversion procedure. E.g., given (5), BESTFIT produces (see Section 3 for details) a ®tted beta distribution ff …x† with ®tted parameters …Uf ; Vf ; pf ; qf † ˆ …ÿ0:451; 23:085; 0:2849; 3:1862† and ®tted €m

…6†

H.-S. Lau et al. / European Journal of Operational Research 120 (2000) 657±670

lf ˆ 1:481; rf ˆ 3:051; b1f ˆ 2:35; b2f ˆ 11:09:

661

…7†

This study evaluates the closeness between ft …x† and ff …x† by considering the di€erences of their €m. Thus, comparing Eqs. (2) and (7), the closeness measures (or ``moment deviations'') are expressed as: Dl ˆ jlf ÿ lj ˆ 0:481; or 48:1% of r; Db1 ˆ jb1f ÿ b1 j ˆ 0:76;

Dr ˆ jrf ÿ rj ˆ 2:051 or 205:1% of r; Db2 ˆ jb2f ÿ b2 j ˆ 5:995:

…8†

If x was not standardized with l ˆ r ˆ 1, then the deviation in rf should be evaluated in terms of ‰…rf ÿ r†=rŠ. Also, lf should then be evaluated in terms of ‰…lf ÿ l†=rŠ; note that the ``obvious'' alternatives of considering either D1 ˆ …lf ÿ l† or D2 ˆ ‰…lf ÿ l†=lŠ are inappropriate, because both D1 and D2 can be made arbitrarily small by using a small r, and D2 by using a large l. Since the bell-shaped beta's (b1 , b2 ) has the small ranges of ``0 6 b1 6 4'' and ``1:8 6 b2 6 9'', deviations in b1f and b2f can be expressed in terms of …bif ÿ bi † regardless of whether x is standardized w.r.t. l and r. The literature on managerial decision models suggests that the €m-deviations presented in Eq. (8) are among the most important criteria for evaluating the closeness between ft (x) and ff (x). For example, in Hillier's (Hillier, 1963) model, after the subjective distributions of the (Ci )s (year-i's cash ¯ow, i ˆ 1 to n) are estimated for a proposed project, if the (li )s and (ri )s of the (Ci )s can be determined, they then enable the l and r (Hillier's primary decision criteria) of the project's NPV and/or IRR to be determined easily. Errors in converting the (Ci )s' elicited fractiles into their true (li )s and (ri )s translate directly into errors in the computed l and r of the NPV or IRR. Similarly, the importance of b1 and b2 of such variables as the (Ci )s and the NPV is explained in, e.g., Arditti (1967) and Kottas et al. (1978). Many other closeness criteria (e.g., a Kolmogorov±Smirnov-type criterion c1 ˆ sup‰Ff …x† ÿ Ft …x†Š) are less directly related to common managerial decision criteria. In some models (e.g., inventory service levels) only closeness at one tail-end is relevant; however, in these situations one should not consider the current problem of ®tting the entire distribution in the ®rst place. 3. Piecewise-linear approximation to the CDF and the BESTFIT software For converting a given set of fractiles to a continuous distribution function, the formation of a Ffÿ1 (p) by simply joining each adjacent pair of plotted fractiles with a straight line (see Fig. 2) is a common (and sometimes the only) procedure described in simulation textbooks (e.g., Law and Kelton, 1982; Banks and Carson, 1984; Thesen and Travis, 1992). This is also the major procedure used by GPSS (Schriber, 1974) to generate non-uniform random variates. A piecewise-linear Ffÿ1 …p† can only be constructed from a set of fractiles that includes x0 and x1 , but the extensive fractile-elicitation literature strongly suggests that subjective estimates of x0 and x1 should not be elicited because they tend to have much larger estimation errors than other fractiles. Nevertheless, in order to evaluate this very well-known conversion procedure, we will assume that x0 and x1 are given (only) when the piecewise-linear procedure is used. That is, in using any given fractile-set fSk g (e.g., fS4 g in (5a), the piecewise-linear procedure will actually be given two more error-free fractiles than the other conversion methods we will be evaluating. Since (as shown below) the piecewise-linear method turns out to be so much less accurate than other methods, this unfair advantage becomes irrelevant. Given any set of fractiles (including x0 and x1 ), the €m of the piecewise-linear Ffÿ1 …p† can be computed with Eqs. (B.4) and (B.2) derived in Appendix B. For example, substituting the fractiles given in (5) (plus x0 ˆ ÿ0.451 and x1 ˆ 23.085) to Eq. (B.4) and then Eq. (B.2) gives the following ®tted €m and hence moment deviations: lf ˆ 1:481; rf ˆ 3:051; b1f ˆ 22:29; b2f ˆ 26:98;

…9†

662

H.-S. Lau et al. / European Journal of Operational Research 120 (2000) 657±670

Fig. 2. A plot of fractiles.

Dl ˆ 0:481; Dr ˆ 2:051; Db1 ˆ 20:73; Db2 ˆ 21:89:

…10†

Recall that the ®tted-beta's parameters (Uf , Vf , pf , qf ) can be easily computed from the ®tted €m given in Eq. (9) using standard formulas (as done in Eq. (3)). The above procedure was repeated for 999 additional test problems, and the average moment deviations were computed for this subsample. The procedure was then repeated for the second subsample of 1000 test distributions to ascertain that the average moment deviations obtained from the two subsamples are the same. The average moment deviations for the aggregate sample (2000 distributions) are then computed for Table 1. This same approach will be used to evaluate the other conversion methods. Since the moment deviations for the piecewise-linear approximation method are huge, we veri®ed the correctness of Eq. (B.4) by: (i) simulation (i.e., generating random variates from Ffÿ1 …p† and computing the variates' sample €m); and (ii) comparing our computed answers with the answers produced by the BESTFIT software. A major cause of the piecewise-linear-method's huge errors is the in¯ated ®tted moments caused by the piecewise-linear-function's positive tail-end (recall that we use only positively skewed distributions here; a negatively skewed distribution will have the same error-magnitudes). We have found that the moment deviations can be substantially reduced by providing erroneous x0 and x1 that constitute a narrower range than the true x0 and x1 . However, the moment deviations after these improvements are still much larger than those of the other conversion procedures, therefore these improvement e€orts are not reported. 3.1. The BESTFIT software For lf and rf , the Palisade Corporation informed us that BESTFIT obtains them using the piecewiselinear model, and we veri®ed this by comparing the (lf )s and (rf )s produced by running the BESTFIT software with those computed by our piecewise-linear computer program. Therefore BESTFIT's average Dl and Dr are the same as those of the piecewise-linear method (see Table 1).

H.-S. Lau et al. / European Journal of Operational Research 120 (2000) 657±670

663

In contrast, for b1f and b2f , it appears that BESTFIT obtains their initial estimates with a linearpiecewise-related procedure, but these initial estimates are improved through a numerical search procedure to obtain the ®nal b1f and b2f . However, we are unable to obtain sucient details from the Palisade Corporation about the procedure to program it, and it is impractical to test the BESTFIT software with 2000 test problems, since the software requires each problem to be keyed-in and solved individually. We therefore ran the BESTFIT software for the ®rst 20 of our 2000 randomly-generated distributions, and the resultant average Db1 and Db2 are given in Table 1. Although they are smaller than piecewise-linear's average Db1 and Db2 , they are still much larger than the average Db1 and Db2 of the other methods. Adding this to their large Dl and Dr (which ARE based on 2000 test distributions), we feel that it is not worthwhile to use a larger sample of test distributions to further evaluate BESTFIT's performance on the problem considered here.

4. Spline ®tting Many (e.g., Winkler, 1967; Hertz and Thomas, 1983) have suggested that a smooth curve can be ``®tted'' to a given set of fractiles. ``Fitting'' here can be interpreted as either: (i) ``joining'' the plotted fractile-points (see Fig. 2) with curve-segments to form an overall smooth curve; or (ii) ®nding the parameter-values that make a single smooth mathematical function pass ``as close as possible'' to the fractiles. Computerized implementation of these two variations are considered in Sections 4 and 5, respectively. Spline ®tting is a standard numerical curve-®tting procedure. For a given set of ``joints'' (the given fractiles in our case), the procedure de®nes one polynomial segment between each successive pair of ``joints,'' and every pair of polynomial segments meeting at a joint satisfy certain speci®ed continuity or smoothness conditions. Although the polynomial segments can theoretically be of any order, ``for historical and other reasons, cubic splines are the most heavily used . . . '' (IMSL, 1987, p. 402). However, while Ffÿ1 …p† should be monotonically non-decreasing, not all cubic spline-®tting procedures will produce a monotonically non-decreasing ®tted inverse cdf for a given set of non-decreasing fractiles. Akima's (Akima, 1970) spline-®tting procedure explicitly attempts to ensure a ®tted curve's non-decreasing monotonicity by minimizing ``oscillations'' (see, e.g., IMSL, 1987, p. 403), therefore, results for Akima's procedure (executed by IMSL's subroutine CSAKM) are reported in Table 1. We have repeated our experiment with other spline-®tting procedures to con®rm that they do not produce smaller moment-deviations than those obtained with Akima's procedure. Computational note: Given a set of (say) n ˆ 7 fractiles (such as in (5), a spline-®tted function is extrapolated to obtain estimates of x0 and x1 . The spline-®tted Ffÿ1 …p† for the enlarged set of (n + 2) fractiles is then used to compute the ®tted €m using Eqs. (B.1) and (B.2). The integration in Eq. (B.1) is done with closed-form spline-integration formulas, which can be derived using the same principle explained in Appendix B.

5. Tocher's curve We now consider ®tting the fractiles to the Tocher (1963) inverse cdf 2

Fsÿ1 …p† ˆ x ˆ a ‡ bp ‡ cp2 ‡ a…1 ÿ p†
ln…p† ‡ bp2
ln…1 ÿ p† with five parameters a; b; c; a b:

…11†

Ideally, a mathematical function used to ®t a given set of fractiles should satisfy the following conditions: (A) it should always be a valid cdf or inverse cdf (clari®ed below); (B) it should be suciently versatile to

664

H.-S. Lau et al. / European Journal of Operational Research 120 (2000) 657±670

cover the required (b1 , b2 )-area; (C) the required ®tting procedure should be simple. The beta cdf satis®es conditions (A) and (B) but not (C), since the beta does not have a closed-form F(x) or F ÿ1 …p† and the required numerical ®tting process would be very tedious. The Weibull (for example) cdf violates (B) (it covers only a (b1 , b2 ) ``line''), but satis®es (A); it somewhat satis®es (C), since the closed-form cdf can be ®tted by a numerical optimization procedure. We consider Fsÿ1 …p† (see Eq. (11)) in this study because: (i) it is a linear function of its parameters and can therefore be ®tted very easily by linear regression (condition C); (ii) it has ®ve parameters and satis®es (B) very well. It satis®es (A) most of the time; i.e., Fsÿ1 …p† gives de®ned and unique x-values for 0 6 p 6 1; however, for some combinations of parameter-values Fsÿ1 …p† may not be monotonically non-decreasing when p is close to 0 or 1. Other distribution functions were also considered; ultimately, presenting Fsÿ1 …p† here is justi®ed by its performance re¯ected in Table 1. One might say that among those considered so far this is the ®rst conversion method that provides reasonable accuracy. After the parameter-values of Fsÿ1 …p† are determined through linear regression, the ®tted €m can be computed with Eqs. (B.1) and (B.2).

6. Linear functions of fractiles that compute the €m Lau et al. (1996) developed ``linear functions of fractiles'' that can estimate very accurately a betadistribution's l and r. E.g., one such lf -computing function is lf ˆ c1 …x0:01 ‡ x0:99 † ‡ c2 …x0:1 ‡ x0:9 † ‡ c3 …x0:25 ‡ x0:75 † ‡ c4 …x0:5 †;

…12a†

where …c1 ; c2 ; c3 ; c4 † ˆ …0:04; 0:11; 0:23; 0:24†:

…12b†

These functions were developed by: (i) randomly generating 2000 beta distributions and determining the fractiles of each generated distribution; (ii) using linear regression to ®nd the ci -values in (say) Eqs. (12a) and (12b) that will minimize the mean square error of the 2000 resultant (lf )s (see Lau et al. (1996) for details). Here we extend Lau et al.'s results as follows: 1. Functions are now developed to compute b1f and b2f as well as lf and rf . 2. Better accuracy is achieved by using unrestricted ci -values. Explanation: Functions such as Eqs. (12a) and (12b) were developed as competitors of other manual lf - and rf -formulas (e.g., the ``PERT'' and the Pearson±Tukey formulas), hence, function complexity was reduced by rounding o€ the ci -values and by requiring each pair of symmetrical fractiles to have the same ci -value. These restrictions necessarily lead to reduced accuracy. In the current problem no manual procedure is envisaged; since a linear function of any complexity will still take negligible computer time, no restriction is placed on the ci -values. For example, applying Lau et al.'s (1996) procedure to our 2000 random beta distributions for the fractileset {S4 } gives the following lf -function with unrestricted ci -values: lf ˆ 0:0784…x0:05 † ‡ 0:1851…x0:10 † ÿ 0:0289…x0:25 † ‡ 0:4100…x0:5 † ‡ 0:3569…x0:75 † ÿ 0:3032…x0:90 † ‡ 0:3017…x0:95 †:

…13†

Appendix A gives the ci -values for the lf -, rf -, b1f - and b2f -functions for all the 10 {Sk }s. Those (ci )s are ``improved'' as explain in Section 7, and di€er from, e.g., the (ci )s illustrated in Eq. (13).

H.-S. Lau et al. / European Journal of Operational Research 120 (2000) 657±670

665

7. Results 7.1. Basic results Computations described above were done for all 10 {Sk }s, their results follow the same pattern depicted by the results for {S4 } and {S3 } in Table 1: the ``Tocher-curve'' (``TC'') and ``linear-fractile-function'' (``LFF'') methods are considerably more accurate than the other three alternatives. Results for all 10 {Sk }s will be presented later in Table 2 in a more relevant situation. Table 1 shows that the CPU time (with an IBM 3090-200S) required for the piecewise-linear method is the second lowest. The CPU time for BESTFIT's numerical optimization procedure beyond the piecewiselinear computations is unknown. However, we feel that their large moment-deviations make the CPUfactor irrelevant. The LFF method is the fastest method, it is also somewhat more accurate than the TC method. Its major disadvantage is that if the given fractiles do not correspond to any one of the 10 {Sk }'s, the (ci )s in Appendix A are invalid, and it is of course impractical then to develop a new set of (ci )s for the given fractiles. However, given the advantages of using one of the {Sk }s to elicit fractiles, it is hoped that given elicited fractiles will often conform with the {Sk }s. Nevertheless, if fractiles at arbitrary fractile-levels are given, the TC method can be used. Although the TC-method's superiority has been proven so far only w.r.t. to the 10 {Sk }-fractile-sets, we have con®rmed its (expected) superiority to the other methods by repeating the computations described in Sections 3±5 using fractiles at randomly generated fractile-levels. TC-method's comparatively high CPU time (see Table 1) is due to Eq. (B.1)'s numerical integration. Recall that for the other conversion methods that require integrations pertaining to Eq. (B.1), the integrations could be performed with closed-form expressions. Nevertheless, our guess is that this CPU time should still be less than that required by BESTFIT's numerical optimization procedure. 7.2. E€ects of a di€erent class of parent distributions Many of the common f(x)s assumed in the management science literature ± e.g., normal, gamma, exponential, Weibull ± have at least one in®nite tail. We indicated in Section 3 that the very common ``piecewiselinear approximation procedure'' requires the speci®cations of x0 and x1 ; thus, this procedure would break down and could not have been tested if in®nite-tailed test distributions were used (because specifying either x0 ˆ ÿ1 or x1 ˆ 1 would give rf ˆ 1 due to the approximating straight-line tail connecting, say, x0:99 to x1 ). Also, each of the above-mentioned common f(x)s corresponds to merely a point or a line in or bordering the beta (Type-1) (b1 , b2 )-area (see Fig. 1). For the above two reasons, we have assumed that the parent distributions are ®nite-tailed beta distributions (as in, e.g., Perry and Greig, 1975; Keefer and Verdini, 1993). However, for a probability distribution with a (b1 , b2 ) in the beta (Type-I) area, the parent distribution may actually follow (say) a Johnson-SB f(x) (see, e.g., Hahn and Shapiro, 1967). We investigate below this Table 2 Average moment deviations for all 10 fractile-sets {Sk }s, from the Tocher-curve (TC) and linear-function (LFF) methods (1000 random beta and 1000 random Johnson-SB test distributions) {S1 }

Dl (% of r) Dr (% of r) Db1

Db2

{S2 }

{S3 }

{S4 }

{S5 }

{S6 }

{S7 }

{S8 }

{S9 }

{S10 }

TC

LFF

TC

LFF

TC

LFF

TC

LFF

TC

LFF

TC

LFF

TC

LFF

TC

LFF

TC

LFF

0.21 0.51 0.09 0.22

0.04 0.37 0.10 0.12

0.10 0.32 0.05 0.17

0.04 0.41 0.11 0.12

0.04 0.20 0.06 0.18

0.04 0.41 0.12 0.11

0.17 0.84 0.12 0.31

0.04 0.41 0.14 0.12

0.37 0.82 0.12 0.26

0.06 0.57 0.13 0.15

0.15 0.29 0.08 0.21

0.04 0.49 0.13 0.16

0.17 0.64 0.07 0.21

0.05 0.44 0.13 0.18

0.23 0.78 0.09 0.26

0.06 0.58 0.13 0.20

0.17 0.79 0.14 0.35

0.06 0.62 0.13 0.23

TC 0.35 1.50 0.19 0.43

LFF 0.11 0.74 0.13 0.26

666

H.-S. Lau et al. / European Journal of Operational Research 120 (2000) 657±670

enlarged problem: A set of fractiles has been elicited from a parent distribution known to have a (b1 , b2 ) in the bell-shaped beta area, but the parent distribution form may be either beta or Johnson-SB ; how will the conversion procedures perform? Out of the earlier 2000 test distributions with given (b1 , b2 )s, only 1000 are now assumed to be beta distributed and hence their fractiles remain unchanged. The other 1000 are now assumed to be Johnson-SB distributed, and their fractiles are now re-computed to correspond to a Johnson-SB distribution. The di€erence between the two parent-distribution assumptions is illustrated by (5c) and (14b). (5c) gives the beta fractiles. With the same €m stated in Eq. (2), if the parent distribution is instead Johnson-SB , the fractiles are: Fractile level ai Fractile value x… i† ˆ xai

0.05 ÿ0.207

0.10 ÿ0.060

0.25 0.267

0.50 0.785

0.75 1.507

0.90 2.353

0.95 2.948

(14a) (14b)

With the modi®ed set of fractiles of 2000 distributions, the procedures described in Sections 4±6 were repeated. Since the spline-®tting method again gives comparatively larger moment deviations, Table 2 presents the moment deviations for all {Sk }s, but only for the TC and LFF methods. Comparing the results presented in Tables 1 and 2 for {S3 } and {S4 } shows that, as expected, the moment deviations deteriorate somewhat when the parent-distribution-type is expanded from ``only beta'' to ``either beta or Johnson-SB .'' However, the extent of deterioration is quite small. Similarly, although {S3 } gives better results than (say) {S10 }, the e€ect of using di€erent {Sk }s is not enormous. The LFF's ci -values in Appendix A are from this latter run and not from the earlier run with pure beta fractiles, therefore {S4 }'s (ci )s in line A of Appendix A are not the same as those in Eq. (13). However, we found that the (ci )s obtained from the pure-beta sample will still perform very well for the ``mixed'' sample (i.e., with Johnson-SB distributions). For distributions in the Type-I (b1 , b2 )-area, considering both the beta and Johnson-SB functions represents a ``reasonable'' coverage of likely distributions typically assumed in the management science literature. However, since a given set of fractiles (as in (5)) can correspond to an in®nite number of parent mathematical distribution functions, it is always possible to ``trick'' any given conversion procedure by devising a f(x) that will lead to large conversion errors. As an extreme example, assume that the (unknown) parent distribution is ``stable Paretian'' ± this distribution was widely advocated for modeling a large variety of ®nancial variables (see, e.g., Fama (1965) and reference list in Lau et al. (1990)), its r ˆ b1 ˆ b2 ˆ 1. Using a right combination of parameters, a stable-Paretian distribution can have ``true'' fractiles that are very similar to (say) those given in (5) or (14), and hence these fractiles would lead most conversion methods to produce some ®nite rf , b1f and b2f , hence Dr ˆ Db1 ˆ Db2 ˆ 1! Thus, it is unrealistic to expect a conversion method to do well for all possible f(x)s. The TC and the LFF methods do not require knowing the parent distribution, neither does the momentdeviation closeness criterion; they all operate via the €m, whose advantage has now been highlighted by the ``mixed sample'' scenario. This scenario becomes very troublesome if either the conversion method or the closeness criterion requires the speci®cation of a parent distribution (e.g., when a Kolmogorov±Smirnovtype criterion is used), because it is often unrealistic to expect the expert providing the subjective fractiles to also specify whether the parent distribution is (say) beta or Johnson-SB .

8. Conclusion and future extensions For converting a given set of fractiles to a distribution function or its moments, we showed that the performance (as measured by ``moment derivations'') of several alternative procedures di€er considerably, and that the most well-known alternative often produces huge errors. The Tocher-curve and the linear-

H.-S. Lau et al. / European Journal of Operational Research 120 (2000) 657±670

667

function-of-fractiles methods are shown to be much more accurate. The €m produced by these methods are often what decision models need. If one chooses to specify a distribution form (e.g., a beta), standard formulas are available to compute the distribution's parameters from the €m. Tocher's Fsÿ1 …p† is particularly convenient for variate generation in simulation. Two meaningful future extensions are discussed below. 8.1. E€ects of errors in the given estimated fractiles The closeness of ff (x) to ft (x) depends on the conversion procedure as well as on the accuracy of the given fractiles. To remove the confounding e€ect of fractile inaccuracy, this study conforms with the earlier related studies (e.g., Pearson and Tukey, 1965; Perry and Greig, 1975; Keefer and Verdini, 1993) and use error-free fractiles. In reality a subjectively estimated fractile contains error that vary with (among other things) the fractile-level and the level of the human-estimator's expertise. It is unknown at this point how these factors should be used to in¯uence the choice of a conversion procedure. 8.2. Using a more general assumption of parent distribution form We indicated earlier that a distribution with a (b1 , b2 ) in the bell-shaped-beta area need not necessary have a beta f(x) as given in Eq. (1). A further point to be noted at this stage is that a real-life distribution's (b1 , b2 ) need not be in the beta area; it can be in the vast area (see Fig. 1) below the ``gamma line'' which forms the lower boundary of the beta area. However, most established management science models do not model random variables with (b1 , b2 )s outside the beta area. By restricting ourselves to that (b1 , b2 )-area assumed by the management science literature as the most relevant, this initial study is then able to evaluate ``fairly'' the most commonly known conversion procedure ± the piecewise-linear method (including the BESTFIT variant). This is because while the Pearson Type-1 (i.e., beta) and Johnson-SB distributions have ®nite tails (i.e., x0 and x1 ), the f(x)s for the (b1 , b2 )-area below the beta area (namely, the Pearson Type-VI and Type-IV and the Johnson-SU ) have at least one in®nite tail; however, the piecewise-linear method can only be implemented after ®nite values are provided for both x0 and x1 . After justifying the exclusion of the piecewise-linear method, the performance of other conversion methods with distributions outside the beta (b1 , b2 )-area can be considered in a subsequent study. Appendix A. Further details of the linear-function-of-fractiles (LFF) method The (ci )s of moment-computing LFFs are given below for each of 10 possible {Sk }s. Lines ``A'' and ``B'' give the (ci )s for the lf - and rf -functions, respectively. Thus, if the {S1 }-fractiles are used, the coecient for rf -functions are ÿ0.56229Eÿ2 and 0.16227E+1, respectively. Line-C coecients x…2† ˆ x0:05 in the lf - and p  6 6 Š, NOT b are for computing x ˆ r‰ b 1f (which is obtained with an additional step of x ). This is because 1f p  6 as l, r and the fractiles. Similarly, the line-D coefb1 is the quantity that has the same dimensionality p ®cients are for a LFF that computes ‰r
4 b2f Š. Our LFFs applies to a positively skewed distribution. If the given fractiles are such that ``…x0:5 ÿ xa † > …x1ÿa ÿ x0:5 †'' for a ˆ 0.01 and/or 0.05 (depending on which {Sk } is used), the parent distribution is negatively skewed. In this case, one obtains an equivalent ``mirror image'' positively skewed distribution with the operation. ‰xa of negatively-skewed f …x†Š ˆ ÿ‰x1ÿa of positively-skewed f …x†Š:

…A:1†

668

H.-S. Lau et al. / European Journal of Operational Research 120 (2000) 657±670

After substituting the converted fractiles from Eq. (A.1) pinto  the appropriate LFFs, switch the negative sign of the computed lf , and remember that now a1f ˆ ÿ b1f . The (Ci )s of moment-computing LFFs

{S1 } ˆ {1 2 3 5 7 8 9} A 0.28114Eÿ1 ÿ0.56229Eÿ2 0.22198E+0 0.50460E+0 B ÿ0.42890E+0 0.16227E+1 ÿ0.16681E+1 0.12570E+0 C 0.13261E+1 ÿ0.53126E+1 0.49237E+1 ÿ0.19569E+1 D ÿ0.10232E+1 0.38509E+1 ÿ0.33474E+1 ÿ0.15971E+0

0.28196E+0 0.51803E+0 0.26327E+1 0.25226E+1

{S2 } ˆ {1 2 4 5 6 8 9} A 0.19894Eÿ1 0.10288E+0 0.26119E+0 0.21554E+0 B ÿ0.23758E+1 0.35151E+0 ÿ0.66454E+0 ÿ0.72564Eÿ2 C 0.15850E+1 ÿ0.38990E+1 0.66842E+1 ÿ0.80227E+1 D ÿ0.31554E+0 0.37321E+0 0.55849E+0 ÿ0.30024E+1

0.29690E+0 0.73533Eÿ1 0.30081Eÿ1 0.48810E+0 ÿ0.10742E+0 0.17710E+0 0.42881E+1 ÿ0.12987E+1 0.66308E+0 0.31913E+1 ÿ0.17916E+1 0.98622E+0

{S3 } ˆ {1 3 4 5 6 7 9} A 0.29699Eÿ1 0.18197E+0 B ÿ0.10409E+0 ÿ0.33265Eÿ1 C 0.95979E+0 ÿ0.51646E+1 D 0.31378Eÿ1 ÿ0.15091E+1

0.96170Eÿ1 0.36020E+0 0.58634Eÿ1 ÿ0.90293E+0 0.93956E+1 ÿ0.94109E+1 0.44727E+1 ÿ0.78118E+1

0.16132E+0 0.13280E+0 0.37860Eÿ1 0.11700E+1 ÿ0.36895E+0 0.18045E+0 0.56851E+1 ÿ0.19551E+1 0.49038E+0 0.76959E+1 ÿ0.37218E+1 0.84210E+0

{S4 } ˆ {2 3 4 5 6 7 8} A 0.13714E+0 0.79036Eÿ2 B 0.94834Eÿ1 ÿ0.10570E+1 C 0.36213E+1 ÿ0.79227E+1 D 0.18392E+1 ÿ0.62941E+1

0.25810E+0 0.90513Eÿ1 0.21482E+1 ÿ0.36779E+1 0.87466E+1 ÿ0.88689E+1 0.11672E+2 ÿ0.17538E+2

0.58891E+0 0.42985E+1 0.77210E+1 0.19751E+2

{S5 } ˆ {1 2 5 8 9} A 0.87133Eÿ2 0.17155E+0 0.63678E+0 0.18319E+0 ÿ0.18103Eÿ3 B 0.18264E+0 ÿ.53430E+0 0.79103Eÿ1 0.19029E+0 0.82305Eÿ1 C 0.25611E+0 ÿ0.24900E+0 ÿ0.32521E+0 ÿ0.27581Eÿ1 0.34590E+0 D 0.46597E+0 ÿ0.94190E+0 0.27922E+0 ÿ0.37157E+0 0.56834E+0 {S6 } ˆ {1 3 5 7 9} A 0.34216Eÿ1 0.19645E+0 0.53678E+0 0.20169E+0 B 0.28128Eÿ1 ÿ0.40276E+0 0.71041Eÿ1 0.18379E+0 C 0.18401E+0 ÿ0.19387E+0 ÿ0.30337E+0 ÿ0.26051Eÿ1 D 0.26483E+0 ÿ0.90424E+0 0.55721E+0 ÿ0.42549E+0

0.30890Eÿ1 0.11982E+0 0.33955E+0 0.50780E+0

{S7 } ˆ {1 4 5 6 9} A 0.52011Eÿ1 0.50366E+0 ÿ0.11437E+0 0.50868E+0 0.50011Eÿ1 B ÿ0.71914Eÿ1 ÿ0.38444E+0 ÿ0.93707Eÿ2 0.32349E+0 0.14220E+0 C 0.13580E+0 ÿ0.22992E+0 ÿ0.18763E+0 ÿ0.53003Eÿ1 0.33508E+0 D 0.12194E+0 ÿ0.17089E+1 0.22577E+1 ÿ0.11380E+1 0.46742E+0

ÿ0.75230Eÿ1 ÿ0.36536E+0 ÿ0.23294E+1 ÿ0.28913E+1

ÿ0.41208E+0 ÿ0.35726E+1 ÿ0.68121E+1 ÿ0.17174E+2

0.44211Eÿ1 0.19597E+0 0.71629E+0 0.10480E+1

0.32952E+0 0.17657E+1 0.35152E+1 0.77434E+1

H.-S. Lau et al. / European Journal of Operational Research 120 (2000) 657±670

{S8 } ˆ {2 3 5 7 8} A 0.19727E+0 B 0.38239E+0 C 0.13398E+1 D 0.21292E+1

ÿ0.16382Eÿ1 ÿ0.91594E+0 ÿ0.17831E+1 ÿ0.34765E+1

0.63387E+0 0.66703Eÿ2 0.40380E+0 ÿ0.55263E+0 0.60138E+0 ÿ0.20628E+1 0.20043E+1 ÿ0.35463E+1

0.17863E+0 0.68252E+0 0.19053E+1 0.28900E+1

{S9 } ˆ {2 4 5 6 8} A 0.18749E+0 ÿ0.33557Eÿ2 B 0.58627Eÿ2 ÿ0.11169E+1 C 0.69351E+0 ÿ0.27194E+1 D 0.84503E+0 ÿ0.50996E+1

0.62292E+0 0.12654Eÿ1 0.14888E+1 ÿ0.86588E+0 0.38740E+1 ÿ0.30102E+1 0.78918E+1 ÿ0.52562E+1

0.18035E+0 0.48835E+0 0.11629E+1 0.16201E+1

{S10 } ˆ {3 4 5 6 7} A 0.51672E+0 B 0.29705E+0 C 0.21335E+1 D 0.26959E+1

0.15152E+1 0.32462E+1 0.87250E+1 0.14406E+2

0.51251E+0 0.13328E+1 0.32159E+1 0.44474E+1

ÿ0.77629E+0 ÿ0.21644E+1 ÿ0.64361E+1 ÿ0.99757E+1

ÿ0.76808E+0 ÿ0.27113E+1 ÿ0.76370E+1 ÿ0.11572E+2

669

Appendix B. Computing moments from an inverse CDF x ˆ F ÿ1 …p† De®ne l0k ˆ E…xk † as the ``kth moment about zero.'' It is known that l0k

Z1 ˆ

‰F ÿ1 …p†Šk dp:

…B:1†

0

Given the (l0k )s, the (lk )s can be computed with the standard formulas (Kendall et al., 1987): l2 ˆ l02 ÿ …l01 †2 ; l3 ˆ l03 ÿ 3l01 l02 ‡ 2…l01 †3 ; l4 ˆ l04 ÿ 4l01 l03 ‡ 6…l01 †2 l02 ÿ 3…l01 †4 :

…B:2†

B.1. Closed-form l0k -formula for a piecewise-linear inverse CDF Fig. 2 shows a plot of fractiles x…k† (k ˆ 1 to n) and their corresponding fractile levels (Pk )s. Each linear segment i (between Piÿ1 and Pi ) of the piecewise-linear function F ÿ1 …p† can be expressed as F ÿ1 …p† ˆ x ˆ x…iÿ1† ‡

x…i† ÿ x…iÿ1†
…p ÿ Piÿ1 † ˆ x…iÿ1† ‡ Si
…p ÿ Piÿ1 †; Piÿ1 6 p 6 Pi ; Pi ÿ Piÿ1

…B:3†

where Si ˆ ‰x…i† ÿ x…iÿ1† Š=‰Pi ÿ Piÿ1 Š is the slope of the ith segment of Fÿ1 (p). Therefore, ZPi ‰F

ÿ1

k

…p†Š dp ˆ

Piÿ1

and l0k ˆ

ZPi

h i k K‡1 ‰x…iÿ1† ‡ Si
…p ÿ Piÿ1 †Š dp ˆ x…i† ÿ xk‡1 …iÿ1†  ‰Si …k ‡ 1†Š

Piÿ1

8P n‡1 < Z i X iˆ1

:

Piÿ1

k

‰F ÿ1 …p†Š dp

9 = ;

ˆ

n‡1 nh X iˆ1

o i k‡1 x…i† ÿ xk‡1 …iÿ1†  ‰Si …k ‡ 1†Š :

…B:4†

670

H.-S. Lau et al. / European Journal of Operational Research 120 (2000) 657±670

References Akima, H., 1970. A new method of interpolation and smooth curve ®tting based on local procedures. Journal of the ACM 17, 589±602. Alpert, M., Rai€a, H., 1982. A progress report on the training of probability assessors. In: Kahneman, D., Slovic P., Tversky A. (Eds.), Judgement under Uncertainty: Heuristics and Biases, Cambridge University Press, New York. Arditti, D., 1967. Risk and the required return on equity. Journal of Finance 22, 42±49. Banks, J., Carson, J., 1984. Discrete-Event System Simulation. Prentice-Hall, Englewood Cli€s, NJ. Fama, E., 1965. Portfolio analysis in a stable paretian market. Management Science 11, 404±419. Hahn, G., Shapiro, S., 1967. Statistical Models in Engineering, Wiley, New York. Hertz, D., 1964. Risk analysis in capital investment. Harvard Business Review 42, 95±106. Hertz, D., Thomas, H., 1983. Risk Analysis and its Applications, Wiley, New York. Hillier, F., 1963. The derivation of probabilistic information for the evaluation of risky investments. Management Science 9, 443±457. IMSL, 1987. IMSL Library User's Manual, IMSL, Houston, TX. Johnson, N., Kotz S., Balakrishnan, N., 1995. Continuous Univariate Distributions, vol. 2. Wiley, New York. Keefer, D., Verdini, W., 1993. Better estimation of PERT activity time parameters. Management Science 39 (9), 1086±1091. Kendall, M., Stuart A., Ord, J., 1987. Kendall's Advanced Theory of Statistics, vol. 1. Oxford University Press, New York. Kottas, J., Lau, A., Lau, H., 1978. A general approach to stochastic management planning model: An overview. The Accounting Review 53, 389±401. Lau, A, Lau, H., Wingender, J., 1990. The distribution of stock returns: New evidence against the stable model. Journal of Business and Economic Statistics 8, 217±223. Lau, A., Lau, H., Zhang, Y., 1996. A simple and logical alternative for making PERT time estimates. IIE Transactions 28, 183±192. Law, A., Kelton, W., 1982. Simulation Modeling and Analysis. McGraw-Hill, New York. Palisade Corporation, 1996. BESTFIT User's Guide, @RISK User's Guide, New®eld, New York. Pearson, E., Tukey, J., 1965. Approximate means and standard deviations based on distances between percentage points of frequency curves. Biometrika 52 (3/4), 533±546. Perry, C., Greig, I., 1975. Estimating the mean and variance of subjective distributions in PERT and decision analysis. Management Science 21, 1477±1480. Schriber, T., 1974. Simulation Using GPSS. Wiley, New York. Selvidge, J., 1980. Assessing the extremes of probability distributions by the fractile method. Decision Sciences 11, 493±502. Solomon, I., 1982. Probability assessment by individual auditors and audit teams: an empirical investigation. Journal of Accounting Research 20, 689±710. Thesen, A., Travis, L., 1992. Simulation for Decision Making. West Publishing, Eagan, MN. Tocher, K., 1963. The Art of Simulation, English Universities Press, London. van Lenthe, J., 1994. Scoring-rule feedforward and the elicitation of subjective probability distributions. Organizational Behavior and Human Decision Processes 59, 188±209. Winkler, R., 1967. The assessment of prior distribution in Bayesian analysis. Journal of the American Statistical Association 62, 776± 800.