- Email: [email protected]
Su models of the maximum of several dependent random variables: Development and computer application
Comput. & Indus. EngngVol. 10, No. 4, pp. 335-346, 1986
0360-8352/86 $3.00 + 0.00 Pergamon Journals Ltd
Printed in Great Britain
Su M O D E L S OF T H E M A X I M U M OF S E V E R A L D E P E N D E N T RANDOM VARIABLES: DEVELOPMENT AND COMPUTER APPLICATION W . E.WILHELM Department of Industrial and Systems Engineering, The Ohio State University, Columbus OH 43210, U.S.A.
(Received 12 June 1985; received for publication 10 October 1985) Abstract - - An approach for modeling the distribution of the random variable (U) that is the
maximum of several dependent random variables is described in this paper. Applications of the model to various problems associated with the Industrial Engineering profession are described and one - - the accumulation of components in small-lot assembly systems - - is studied in detail. Numerical tests investigate the accuracy of the approach and identify several fundamental characteristics of the accumulation process. Apparent advantages of the model are that it requires computer run time that is a fraction of that required by Monte Carlo sampling (simulation) procedures and that it is readily amenable to parallel processing.
1. INTRODUCTION
A number of models in the field of Industrial Engineering require definition of the probability density function of the random variable that is the maximum of several dependent random variables. PERT networks, for example, require such a maximum to define the "realization" time of a node, the time at which all immediately preceding tasks are completed. According to Whitehouse[1], GERT network analysis entails the maximum operator to model AND nodes. More recently, Wilhelm and Ahmadi-Marandi[2] described an approach, based on the maximum operator, for modeling the transient performance of production/assembly queuing networks. Operation finish times are assumed to be related by the multivariate normal distribution and operation start times are then defined by evaluating the maximum of several dependent, normal random variables: Operation Start Time = max[X1, X2, X3] in which X l = finish time at the upstream station that readies the part, X2 = finish time at the upstream station that readies the subassembly, X 3 -- finish time of the operation that readies the assembly station. A similar evaluation has been used to estimate lifetimes of reliability systems in which components exhibit normally distributed lifetimes: System Lifetime = min [Xl, X2] = - m a x [ - X l , -X2] in which Xl and X2 are lifetimes of components 1 and 2 (resp.). Motivated by these applications, the objective of this study was to investigate the computational advantages of the S, distribution as a model of the random variable that is the maximum of several dependent random variables. The Su distribution, developed by Johnson[3], was selected for 2 reasons. First, the distribution is related to both the normal and lognormal distributions, since it is based on a transformation of a normally distributed random variable. Secondly, the distribution can be fit to a broad range of cases; in fact, it can fit all ([31, [52) combinations (i.e. measures of skewness and kurtosis) which lie in the half-plane below the "line" of ([51, [32) values that can be fit by the ~A~ ~0:4-E
335
336
W.E. WIHIEI,M
lognormal distribution. Based on the theory of extreme value distributions (e.g. see Johnson and Kotz[4]) it is expected that this region is the one of primary interest in practical applications. Prior theoretical research by mathematical statisticians (which is summarized by David[5], Galambos[6] and Gumbel[7]) has led to description of families of extreme value distributions, the limiting distributions (as J---> 2) of the maximum of J i.i.d. random variables. Most of the related, applied research has been motivated by the development of PERT and is summarized by Elmaghraby[8]. A variety of modeling approaches have been proposed including polynomial approximations by Martin[9], bounding distributions by Kleindorfer[10] and Shogan[ll], and hybrid schemes combining numerical and sampling techniques by Hartley and Wortham[12] and Ringer[13]. Clum[14] has recently investigated the use of Burr and S-D distributions to analyse PERT networks. None of these approaches offers the combination of accuracy and computational speed required in the applications described above. Clark[15] developed an alternate approach using the normal distribution as a modeling convention. He derived the (exact) first 4 moments of the random variable which is the maximum of 2 random variables that are related by the bivariate normal distribution. His convention is based on the assumption that this maximum is normally distributed so that the maximum of J dependent random variables can be modeled using a ( J - 1) stage (i.e. pairwise calculation indicated by max[Xl, X2 . . . . . XI] = m a x [ X j _ l . . . . max(X3, max(X2, Xl) • • .)]. This convention was applied to model assembly systems by Wilhelm and Ahmadi-Marandi[2] and was found to give "reasonable" estimates of system performance. Nevertheless, it is expected that accuracy can be improved by using conventions that permit more faithful representation of the maximum by a positively skewed distribution that more closely approximates third and fourth moments than does the (symmetrical) normal distribution. Wilhelm[16] developed and tested a convention based on the lognormai distribution, which is, of course, skewed and was shown to give accurate results with computation time a fraction of that required by sampling procedures. Results of the current research extend this line of investigation to define models based on the S, distribution. It is expected that this new convention will provide a fitting capability important in many applications. The primary advantage of these conventions is that they facilitate computer modeling of events that are defined by the maximum operator. Simulation can be applied to model these events, but runtime may be excessive and results may require involved statistical interpretation. Run times to implement the conventions are typically a fraction of the times that equivalent simulation models require; this attribute is particularly important for effective use of personal computers. In addition, the conventions can be implemented using parallel processors much more readily than simulation approaches. Finally, conventions offer convenient conceptual models which can lead to analytical models to replace simulation as the only available means of modeling certain processes. The body of this paper is organized in 3 sections. The S, convention, presented in Section 2, may be used to approximate the moments of the maximum of multivariate S, random variables, max[Xi, X2 . . . . . Xs], by taking the variables pairwise: U l = X 1
Ui = max[Xj, Uj-l] forj = 2, 3 . . . . .
J.
Equations which state the moments of Uj are given in Section 2 along with the product moment E[X3, U2] which is used to calculate the correlation r[Xj, Uj-*] required to estimate the moments of U. In addition, a procedure for fitting parameters to U~ is given. The convention is demonstrated in Section 3 in an application which involves coordinating the accumulation of components for assembly. This is an important practical problem which determines the amount of in-process inventory as well as schedule (due
S. models of the maximum of several dependent random variables
337
date) performance that can be achieved by an assembly system. Conclusions are given in Section 4. It should be noted that all derivations are outlined in the Appendices. 2. THE S, CONVENTION The which Z ~,, 6 Y
S, distribution is defined, according to Johnson [3], by Z = 5' + 6 s i n h - l ( Y ) , in = a unit normal variate = parameters of the S, distribution = the random variable described by the S, distribution
Equivalently, Z may be translated (i.e. N = ~t + o Z ) to yield the normally distributed random variable N with mean ~t and variance o 2, which, in turn, may be transformed to obtain Y, Y = sinh (N) = (e N - e-N)/2; the corresponding inverse transformation is N = sinh-l(Y) = f n [ Y + (y2 + 1)1/2]. The random variable Y that is defined by the transformation is said to have the 2-parameter S, distribution. The associated normal distribution, related by the transformation sinh(N), has parameters ~t = -5'/6 and o 2 = 1/6 z. If YI and Y2 have the bivariate S, distribution with coefficient of correlation Q(Y1, Y2), the correlation between the associated normal random variables Pl2 can be shown (see Appendix A) to be: Q,z = (1/olo2) tn{[S1 + (Sl 2 + 4S2S3)1/2]/(2S2)}
(1)
in which V[-] and E[.] are used as the variance and expected value operators (resp.) and S, = {o(Y,,
Y2)v~/z[Y,IVm[Y2I + E[Y, IE[Y2I)4
e x p [ - ( o l 2 + 022)/2]
Sz = exp(tx2 + Ix1) + exp(-Ix2 - gl) $3 = e x p ( - g 2 + gl) + exp(g2 - ~q). Flexibility may be enhanced by using the 4-parameter S, distribution to allow fitting the first 4 moments of a distribution, Based on a translation of Y, the random variable X = kY + ~ has mean E[X] = kE[Y] + ~ and central moments vi(X) = kivi(Y) for i >/2. It is easy to show that X and Y h a v e the same ([31, [32) values and that o(Yl, Y2) = 9(X1, X2). The 2-parameter case results from setting ~l = ~ 2 = 0 and kt = ~.2 = 1. A convention based on the S, distribution offers an appealing flexibility, since parameters of this distribution can be set to cover all points in the ([31, [3z) plane that lie "below" the lognormal (L) "line" (Johnson and Kotz [17]). For y < 0 the S, distribution is skewed to the right and may, therefore, provide a good fit for the distribution of max(Xl, )(2, . • . , X:) in which the Xj{j I J = 1, 2 . . . . . J} have the multivariate S,, distribution.
Moments of the maximum of bivariate Su random variables The i tn moment about zero for U = max(X1, X2), Ui, may be calculated for the special case (see Appendix B) in which El = ~ 2 and ~tk > ln(8kk'/~.k) + ok Z~ (for k = 1 and 2). Zk in the latter restriction may be replaced by 5, since Pr[ I Zk I <5] is close to 1.0 (for k = 1 and 2). These restrictions tend to limit the inherent flexibility of the 4-parameter distribution but may be satisfied in a number of situations of practical interest. In applications for which the restriction is reasonable, {/)i ] i = 1, 2, 3, 4} are defined (see Appendix B) by 2
0:
}-" [k,I/Clk + ~,F,] k:l 2
02 = 2 k=l
[kkW2k + 2~kkkl'9~k + ~2fk]
338
W.E. WILHELM 2
[kkW3k + 3~k;~kW2k + 3Xk ~kWlk
. =
+
~3fk ]
k--I 2
[)"kW4k + 4~k~'kW~k + 6~k~'kW2e +
=
kfk]" (2)
k=l
Terms in eqn (2) are defined in the following paragraphs. First, f,
=
-
.z
+
en(XdX2)l/a}
in which = distribution function of the standard normal a = [012 "Jr- 022 -
2~12G10"2] 1/2.
Due to symmetry, f2 =
in which the operator Q[.] is used to indicate that f2 can be obtained from fl by changing all subscripts 1(2) to 2(1) and will be used to denote a number of symmetrical relationships throughout this paper. If, in fact, ~ = ~2, 2 k=l
and the expressions for/..-/i in eqn (2) can be simplified appropriately. In order to define the l,~'ik terms in eqn (2), it is necessary to compute 4 additional terms, 3;/n, a;/il, a;/i2 and a5/i2: /~'/il = /-7'1 (I){[~1 -- ~'2 q-
iol (ol
-
Q12(J2) q'-
fn(~q/~.2)]/a},
(3)
in which /7,l = exp(ibtl + i2o12/2 ),
(4)
a;/il is/1;/n with i replaced by - i in both eqns (3) and (4), a;/i2 = f2[Ma], a;/i2 = fa[a;/i,] and L2 = f2[£.]. Using these terms, l'~'i~ (i.e. 14'ik for k = 1) is
l,~'i, = 2-i[1(4i, + (-1)iMi, + t,,] for i = 1, 2, 3, 4
(5)
in which
tll = O, t21 = - 2 , tBz = -614/11 and tal = -161'f'21 - 2; and, due to symmetry, 14'~2 = f~[W,I and ti2 = Q[t~a] for i = l, 2, 3, 4. It is interesting to note several relationships implicit in eqns (2)-(5). First, in the 2 parameter case [i.e. with 21 = 22 = 1.0 in eqn (3)] the ith m o m e n t about zero, if'i, of W = max((Y1, Y2) is
l'f'i = l'Vi, + I~i2.
(6)
Secondly, the ith moments about zero, M i, of M = max(L1, L2), in which L~ and L 2 a r e bivariate lognormals, is (see Wilhelm[16]) m i ~-- J~/il -{'- J~/i2
(7)
[still using ~-I = ~,2 = 1.0 in eqn (3)]. Finally, £1 (£2) in eqn (4) defines the i th moments about zero for £ | (£2), assuming 0l = 02 = 0.0 (see Johnson and Kotz[4].
S, models of the maximumof several dependent random variables
339
To summarize, the moments of U = max(X1, )(2) may be determined by computing, in order, (i) (ii) (iii) (iv) (v)
P12 using eqn (1), £1 and/]2 using eqn (4), Mil, ~/i2, /~/il a n d / ~ / i 2 using eqn (3), I~il and ff'n using eqn (5), Ui (for i = 1, 2, 3, 4) using eqn (2).
Product moment Specification of the Su convention is completed by defining the product moment E[X3-max(Xl, X2)] = E[X3. U] (see Appendix C): E[X3.max(X1, X2)] = Tit + Ti2
(for i = 1 only)
in which Til = ~-I~3V[/il jr_ ~l~3fl
+ (~.,K3)(D/4)[E,(1,1) + E l ( l , - 1 ) + E l ( - 1 , 1 ) + E , ( - 1 , - 1 ) ] + (~,k3)(O/2)[/~l(1,1) + / ~ ( 1 , - 1 ) + /~,(-1,1) + E , ( - 1 , - 1 ) ]
(8)
Ti2 = •[Til] and terms are defined in the following paragraphs. First, certain functions of the correlations are defined as R = 1 - P122 D = [(1 - ~122 - QI32 - Q232 q- 2QlaOI3Q23)/2R] G1 = (Q13 - ~12Q23)/R B12 = ol -- Q)1202 and G2 = ffaiG,], B2, = ffa[B,2]. The basic function Ek(',') is:
E , ( K , L ) = K L exp[IK~t3 + Lbt, + (P,a/2a 2) + QI] @[(~q - ~t2 + P,)/al in which, f o r K =
+landL
= _+1,
I=1 Al = IL01 + Ko3(Gt + (112G2) P1 = A1BI2 - GeKR0203 + ~n(~.i/~.2) Q1 = exp{(R/2a2)[Al2022 + BI22G22032 + 2A1Bt2G2K°2°3] en(?,,/?,2)[2P, - ln(~.,/~.2)]/R} -
By symmetry, A2 = Q[A,], P2 = ~2[Pl], Q2 = ~[Q,] and E2(K,L) = ffa[E,(K,L)]. It is of interest to note that Ee defines the product moment in the 2 parameter S,, case (with ~.l = ~.2 = 1.0): E[Ys.max(Yl, I12)] = E[Y3.W] = (D/4)[E,(1,1) + E , ( 1 , - 1 ) + E , ( - 1 , 1 ) + E , ( - 1 , - 1 ) + E2(1,1) + E2(1,-1) + E2(-1,1) + E 2 ( - 1 , - 1 ) ] .
(9)
Given these definitions,
E1(K,L) is E I ( K , L ) with I = 0 a n d / ) 2 ( K , L ) = f2[/~,(K,L)]. Parameter fitting procedure Implementation of the Su convention entails calculating the moments of U and then fitting the parameters of an S,, distribution to these moments. S, fitting procedures have
340
w.E. WII.HECM
been developed by Draper [18], Leslie [19], Johnson and Kotz [4], and more recently by Wilson [20, 21], but none of these is useful in the current application. A new procedure is, therefore, described in this subsection. If ~ and k are fixed (e.g. as they are in the 2 parameter case), the parameters (~t, o) of U may be determined by equating the moments of U [/_71 and U2m defined in eqn (2)] to those of an S. distribution: E[X ~] = /)~ in which,
E[X] = (k/2)e °-~/2[e ~' - e - " l + ~
(10)
E[X 2] = 0~2/4)[e 4°2/2 (e 2~t -- e -2~t) -- 2] + ~2 + ~.eo2/Z[e~t _ e - ' ] .
(11)
and
Equations (10) and (11) may be solved simultaneously for Su parameters ~t = In(F) and o 2 = 2 ln(r)
(12)
in which
~"
=
F 1/2, F = [~'I(F-1)][2(/21 - ~)/k] 2
and F is determined as a root of the polynomial equation F4-vl F3+6F 2-rlF+
1 =0
(13)
in which Vl = 4 + [2(/)1 - ~)/k]4[(4/)2/k 2) - (4~/L2)(2/)1 - ~) +
2]-'.
In the 4 parameters case ~ is fixed by the requirement that ~ = El = ~2, but parameters (k, tl, o) may be fitted to the first 3 moments of U, {/)i [ i = 1, 2, 3} defined by eqn (2). The third moment of X is E[X 3] = (~.3/4)e9°~-/2sinh(3g) + (~2 _ k2/4)(3?~eO-V2)sinh(g)
+ (3~k2/2)e 4''-/2 cosh(2~t) + (~2 _ 3~2/2)~,
(14)
so that (k, tx, o) may be estimated by solving eqns (10), (11), and (14) simultaneously. This approach requires solution of a high order polynomial and is, therefore, replaced by the algorithm stated below. (vi) S e t ~ = 1, k t = k m , B ~ = ~tm, Oe = o m , e > O . (vii) D ~= dE[X3]/dk = (0.75)(~f)2exp[4.5(ot)2]sinh[3~t ~] + (1.5)(k~)[3~ - (2 + /)1)]. = (viii) k t + l = k ~ + Ak. (ix) (x) Use ~ and U +1 in eqn. (12) to calculate ~V+1 and C +~. (xi) Use ~, U +1, I,d +1 and C +l in eqn (14) to calculate E[X3] t+l. (xii) If (1 - e)/)3 -< E[X3] c+1 < 1)3 go to step (viii). Else, £ <-- ¢J + 1 and go to step (vii). (xiii) Fitted parameters are e, )~+l, f + t , C + t . Calculate {E[Xi] I i = 1, 2, 3} and stop. It is implicitly assumed that the parameters of U are not dramatically different from those of Xm in which m is the index (i.e. 1 or 2) of max(~tl, ~t2). To begin, the parameters of U are set equal to those of variable m. Equation (14) is then used to improve the current estimate of ~., U. The rate of change of E[X 3] as a function of ~. (using current values ~tt, or), D t, is computed at step (vii) and used in step (viii) to estimate, according to a first-order Taylor series expansion, the increment Ak that would set E [ X 3] = 03 as desired. ~t' is updated using this increment [step (ix)] and then used, along with ~, to update estimates of (~t, o) using eqn (12). Revised parameters are used in eqn (14) to compute the corresponding value of/~[X3]. If this value is within E of U3, the procedure
S. models of the maximum of several dependent random variables
341
terminates; otherwise, it will iteratively improve parameter estimates until E[X1] = 0~, E[X2] = 02, and E[X3] is acceptably close to/-)3. 3. EXAMPLE APPLICATION
Components required to assemble a product must be accumulated before assembly can begin. For example, a printer circuit board cannot be assembled until the board and all electronic components have been accumulated; i.e. "kitted". Some products such as automobiles may require a kit for each subassembly as well as an accumulation of subassemblies before initiating final assembly. It is well known that the bulk of manufacturing is currently done in small lots to satisfy individual customer preferences. In small-lot assembly systems, components may be specially ordered rather than inventoried, because requirements are unique and products themselves are special ordered. Coordinating the input of materials from vendors and from production facilities is thus an important material management problem in small-lot assembly systems. By starting with the due date of a lot of assembled products and "backing off" an appropriate duration, the time at which an accumulation will be needed to initiate assembly can be readily determined. Material managers must then set a due date for the delivery of each component from its supplier. Due dates are thus used to coordinate material flows from vendors and production facilities to assembly operations. Since delivery time is a random variable in actual systems, component due dates must be set to hedge against the dual risks of early and late delivery. If a component is delivered well ahead of the time at which assembly is to begin, schedule performance will be good, but in-process inventory will be high. Conversely, if a component is delivered late, in-process inventory will be low, but schedule performance will be poor. Due dates must, therefore, be set to provide an appropriate balance between the risks of early and late delivery. The S, convention may be applied to study the accumulation process by using the model, with U ~ = X~ and Uj = max[U j - l , Xj] f o r j = 2, 3 . . . .
,J
in which X~ = delivery time of component j (j = 1, 2 . . . . . J) Uj = time at which component accumulation is complete. Computation procedures given previously in this paper may be used in the following manner to estimate the moments of accumulation time: (xiv) Set j = 1; (xv) Use procedure (i)-(v) to calculate moments of/_fl; (xvi) Use eqn (8) and standard indentities to calculate correlations: r[Xk, Uj] = r[Xk, max(U j - l , Xj] for k = j + l , . . . , J; (xvii) Use procedure (vi)-(xiii) to fit parameters of the S, distribution to moments calculated in step (xv) above. (xviii) I f j = J stop, moments of Uj have been calculated. Otherwise, s e t j = j + 1 and go to step (xv). Systematic use of step (xvi) yields pairwise correlations that are subsquently required in step (xv) to calculate moments. Applying the model in the context of component accumulation, E[Xi] is interpreted as the due date for component j and V[Xj] as the variance of delivery time. This model can be used to gain fundamental insights into characteristics of the accumulation process and could be extended, for example, to develop optimal due dates (although such work is beyond the scope of this paper).
342
W.E.
WII.ItELM
Pairwise correlations among delivery times (i.e. r[Xj, X / , ] ) are zero if the components are produced by independent vendors. However, if they are produced in a flexible manufacturing system (or some other common production facility) pairwise correlations will be nonzero. In an actual application, correlations could be estimated from historical experience or by modeling (e.g. see Wilhelm and Ahmadi-Marandi [2]). In this paper, 3 cases are considered for demonstration p u r p o s e s - Q(X1, X2) = 0.0, 0.15 and 0.30. The marginal mean and variance of component delivery times were assumed to equal 100 (~j = 0.0, ~.j -- 1.0, 5'j = -53.0683, and 6j = 10.0254) in all test cases reported in this section. These values were selected purposely, since they define a "worst case" in which it is most difficult to estimate the moments of Uj. If one delivery time had a mean that is much larger than other means, that delivery time would dominate the definition of accumulation time and make it easier to develop good estimates. Practical experience indicates that it is more difficult to meet assembly schedules if the number of components required in the accumulation J is large. Characteristics of the accumulation process are therefore investigated for 5 levels of J: 2, 3, 4, 5 and 10. In all, 15 cases were tested (3 levels of correlation for each of the 5 levels of J) by both the 2 and (restricted) 4 parameter S , distributions. The S, convention is based on the assumption that the Xj are related by the multivariate Su distribution; it imposes the approximation that Ui are also multivariate S, distributed. This introduces a potential source of error, because the true distribution of the maximum is not known. Since the S,, distribution can, in principle, be fit to a broad range of moments, the error may be of acceptable magnitude. Accuracy of the approach was, therefore, studied by comparing test results with estimates developed using a simulation model. Analytical models of the accumulation process are not available. Comparisons with the simulation model were made using several measures: -
-
KS -- KS test statistic related to the null hypothesis that the distribution function
developed from the simulation was that of the S, distribution (critical value for 30 cells and a = 0.50 is 0.241). - - A = (simulated average of E[UJ]) - (recursion estimate of E [ U ] ) . t = t-test statistic related to H0: A = 0, HA: A 4: O. _ _ ~2 = ~2_test statistic related to Ho: (simulated V[UJ]) = (recursion estimate of V[UJ] HA: (simulated V[UJ]) 4: (recursion estimate of V[UJ]). - - Computer runtimes (all tests were made on a PDP 11/44). -
-
Table 1. S,, convention
0
Case J
Average
Simulation* SD CPU (sec)
~t = ~2 = 0.0
Mean
~,1 = ~,2 = 1,0
Approximation SD CPU (sec)
KS
Results A
0.00
2 3 4 5 10
105.582 108.443 110.626 112.059 116.166
8.61 7.86 7.48 7.10 6.43
142.22 249.95 272.78 445.62 422.92
11)5.623 1118.578 110.549 112.001 116.155
8.74 8.04 7.58 7.24 6.29
0.150 0.285 0.751 1.500 2.784
0.014 0.028 0.018 0.022 0.031
-0.041 -0.142 0.077 0.058 11.1)11
0.15
2 3 4 5 10
105.125 107.712 109.696 111.073 114.859
.8.85 8.31 7.90 7.67 7.31
250.73 215.60 302.68 364.96 665.50
105.183 107.901 109.714 111.063 115.010
8.99 8.46 8,11 7.86 7.11
0.651 0.299 0.582 1.481 4.868
0.020 0.023 0.023 0.036 0.031
-0.1158 -0.190 -11.018 0.010 -0.151
0.30
2 3 4 5 10
104.619 106.923 108.696 110.002 113.482
9.06 8.75 8.32 8.25 8.16
151.21 190.23 264.37 384.02 555.26
104.71)2 107.157 108.828 110.091 113.850
9.23 8.84 8.58 8.38 7.78
0.151 0.351 0.466 0.683 2.800
0.016 0.022 0.025 0.033 0.032
-0.083 -I).234 -11.132 -0.089 -0.368
* 1000 replications. A = (Simulation average) - (approximation mean). KS = KS-test statistics based on 30 cells, S, distribution.
S,, m o d e l s o f t h e m a x i m u m o f s e v e r a l d e p e n d e n t r a n d o m v a r i a b l e s
343
Table 2 . 4 parameter S, convention (~l = ~2) Case J
Average
0.00
2 3 4 5 l0
105.582 108.443 110.626 112.059 116.166
8.61 7.86 7.48 7.10 6.43
142.22 249.95 272.78 445.62 422.92
105.623 108.585 110.549 112.001 116.153
8.73 8.(}4 7.58 7.24 6.28
0.265 0.649 0.817 1.683 3.831
0.014 0.028 0.018 0.022 0.029
-0.041 -0.142 (}.077 0.058 0.013
0.15
2 3 4 5 10
105.125 107.712 109.696 111.073 114.859
8.85 8.31 7.90 7.67 7.31
250.73 215.60 302.68 364.96 665.50
105.183 107.902 109.714 111.063 115.008
8.99 8.46 8.11 7.86 7.11
0.217 0.649 0.548 0.833 6.648
0.020 0.023 0.023 0.036 0.031
-0.058 -I}.190 -(}.1}18 0.010 -0.149
(}.30
2 3 4 5 10
104.619 106.923 11}8.696 I10.(X}2 113.482
9.06 8.75 8.32 8.25 8.16
151.21 190.23 264.37 384.02 555.26
104.702 107.157 108.827 110.090 113.848
9.23 8.84 8.58 8.38 7.78
0.151 0.351 0.435 0.583 3.850
0.016 0.022 0.025 0.033 0.032
-0.083 -0.234 -(}.131 -0.088 -0.366
0
Simulation* SD CPU (see)
Mean
Approximation SD CPU (sec)
KS
Results A
* 1000 replications. A = (Simulation average) - (approximation mean).
KS = KS-test statistics based on 30 cells, S. distribution.
Results for the 2 and (restricted) 4 parameter models are given in Tables 1 and 2 respectively. In all cases, the fitting procedure (step xvii) converged in one iteration (with = 0.001). None of the KS or t-statistics are significant (~ = 0.05) indicating that accumulation time, UJ, is approximately St, distributed and that the estimate of mean accumulation time, E[UJ], is acceptably close to the estimate devloped by the simulation model. All A values are small in absolute value, again indicating the accuracy with which the mean of accumulation time is estimated. X2 statistics are significant at an c~ = 0.05 level only for the variances V[UJ] in the 3 cases in which Q(X1, X2) = 0.30 and J = 10. Variance of accumulation time is thus estimated accurately in all but the most extreme cases tested. In all cases, runtime for the approximation model was a small fraction of that required for the simulation model to complete 1000 replications. The approximation model required approximately the time required to complete 5-7 simulation replications. In addition to these statements concerning the accuracy of the approximating procedure, it is possible to interpret results in relationship to the accumulation process. First, it is noted that a deterministic estimate which does not incorporate the variance of delivery time V[XI], significantly underestimates the time at which an accumulation will be complete. The deterministic estimate would be 100 in all cases and the true mean time of accumulation E[U] is seen to increase as the number of components J increases. Conversely, the variance of accumulation time decreases with an increase in J. This indicates that accumulations will be completed consistently at later times as J increases. Safety lead time, the difference between scheduled assembly time and component due date must increase as J increases to maintain an acceptable balance between the risks of early and late delivery. Results for different levels of pairwise correlations between component delivery times indicates that mean accumulation time decreases (and the variance of accumulation time increases) with an increase in pairwise correlation. This result suggests that assembly schedules may be more easily attained if components are produced in a common facility such as a flexible manufacturing system. This relationship has, apparently, not been reported previously in the literature and may suggest further investigation to understand advantages of integrated production/assembly facilities. 4. C O N C L U S I O N S
An approach based on the S,, distribution is presented in this paper as a means of modeling (approximately) the random variable that is the maximum of several dependCAIE
10:4-F
344
w . E . WILttEI.M
ent r a n d o m v a r i a b l e s . A p p l i c a t i o n s in which this a p p r o a c h might be used include m o d e l i n g t r a n s i e n t q u e u i n g n e t w o r k s (e.g. p r o d u c t i o n / a s s e m b l y s y s t e m s ) , P E R T a n d G E R T n e t w o r k s , a n d s y s t e m reliability. A p p l i c a t i o n to m o d e l the a c c u m u l a t i o n of c o m p o n e n t s in s m a l l - l o t a s s e m b l y s y s t e m s is s t u d i e d in detail. T h e a c c u m u l a t i o n p r o c e s s is i m p o r t a n t to m a t e r i a l flow m a n a g e m e n t , since it d e t e r m i n e s i n - p r o c e s s i n v e n t o r y levels as well as a s s e m b l y s c h e d u l e p e r f o r m a n c e . N u m e r i c a l tests i n d i c a t e t h a t the f u n d a m e n t a l a s s u m p t i o n u n d e r l y i n g the a p p r o a c h - that the m a x i m u m of m u l t i v a r i a t e S , r a n d o m v a r i a b l e s has the S , d i s t r i b u t i o n - - gives a g o o d a p p r o x i m a t i o n to the true d i s t r i b u t i o n . F u r t h e r m o r e , the a p p r o x i m a t i o n p r o c e d u r e gave e s t i m a t e s o f the m e a n of the m a x i m u m of d e p e n d e n t r a n d o m v a r i a b l e s which c o m p a r e f a v o r a b l y with t h o s e d e r i v e d f r o m a s i m u l a t i o n m o d e l . T h e s e results indicate that the a p p r o a c h gives a v i a b l e m o d e l for use in the a p p l i c a t i o n s d e s c r i b e d a b o v e . Test results also l e n d insight into the c h a r a c t e r i s t i c s of the a c c u m u l a t i o n p r o c e s s , i n d i c a t i n g that d e t e r m i n i s t i c p r o c e d u r e s for setting d u e d a t e s are i n a d e q u a t e a n d defining i m p o r t a n t r e l a t i o n s h i p s b e t w e e n the n u m b e r o f c o m p o n e n t s to be k i t t e d a n d the t i m e at which a c c u m u l a t i o n can be c o m p l e t e d . R e s u l t s also i n d i c a t e that a s s e m b l y s c h e d u l e s can be a t t a i n e d m o r e c o n s i s t e n t l y if c o m p o n e n t s are p r o d u c e d in a a c o m m o n facility. This p a p e r suggests f u r t h e r r e s e a r c h to d e v e l o p a l t e r n a t e a p p r o x i m a t i o n s for U as d e f i n e d in e q n (18) to relax r e s t r i c t i o n s which define a c c e p t a b l e v a l u e s o f rt a n d o. A d d i t i o n a l l y , e x t e n d e d m o d e l s of the a c c u m u l a t i o n p r o c e s s m a y be d e v e l o p e d to d e t e r m i n e o p t i m a l d u e d a t e s to b e t t e r c o o r d i n a t e m a t e r i a l flow in small-lot a s s e m b l y systems. P e r h a p s the m o s t significant c o n c l u s i o n s r e l a t e to c o m p u t e r a p p l i c a t i o n s . A l t h o u g h m a t h e m a t i c a l l y i n v o l v e d , the a p p r o a c h gives g o o d e s t i m a t e s of p e r f o r m a n c e with a r e l a t i v e l y low c o m p u t e r r u n t i m e . It is e x p e c t e d that this a d v a n t a g e will p r o m o t e use of s o m e m o d e l s o n p e r s o n a l c o m p u t e r s a n d s o l u t i o n of large-scale p r o b l e m s on larger c o m p u t e r s . G i v e n p a r a m e t e r v a l u e s , a n u m b e r o f e q u a t i o n s can be s o l v e d s i m u l t a n e o u s l y [e.g. steps (xv) a n d (xvi)] so the a p p r o a c h is well s u i t e d for p a r a l l e l processing. Simulation is t h e o n l y p r a c t i c a l a l t e r n a t i v e m e t h o d in m a n y a p p l i c a t i o n s , but it m a y r e q u i r e excessive r u n t i m e a n d is n o t r e a d i l y a d a p t e d to p a r a l l e l processing. Acknowledgements - - This material is based on work supported by the National Science Foundation under
grant no. MEA-8213196. The author gratefully acknowledges the work of Mr Surendra Saboo whose careful programming effort and thoughtful analysis of implementation requirements led to the numerical results reported in this paper. Also, constructive comments of the referees led to improvements in an earlier version of this paper and are hereby acknowledged.
REFERENCES 1. G. E. Whitehouse. Systems Analysis and Design Using Network Techniques. Prentice-Hall, Englewood Cliffs, NJ (1973). 2. W. E. Wilhelm and S. Ahmadi-Marandi. A methodology to describe operating characteristics of assembly systems, liE Trans. 3(14), 204-213 (1982). 3. N. L. Johnson. Systems on frequency curves generated by methods of translation. Biometrika 36, 149-176 (1949). 4. N. L. Johnson and S. Kotz. Distribution in Statistics: Univariate distributions - - 1. Wiley, New York (1970). 5. H. A. David. Order Statistics. Wiley, New York (1981). 6. J. Galambos. The Asymptotic Theory o f Extreme Order Statistics. Wiley, New York (1978). 7. E. J. Gumbel. Statistics o f Extremes. Columbia Univ. Press, New York (1958). 8. S. E. Elmaghraby. Activity Networks. Wiley, New York (1977). 9. J. J. Martin. Distribution of time through a directed acyclic network. Ops Res. 13(1), 46-66 (1965). 10. G. B. Kleindorfer. Bounding distributions for stochastic acyclic networks. Ops Res. 19, 1586-1601 (1971). 11. A. W. Shogan. Bounding distributions for a stochastic PERT network. Networks 7,359-381 (1977). 12. H. O. Hartley and A. W. Wortham. A statistical theory for PERT critical path analysis. Mgmt. Sci. 12(10), B469-B481 (1966). 13. L. J. Ringer. A statistical theory for PERT in which completion times of activities are inter-dependent. Mgmt. Sci. 17(11), 717-723 (1971). 14. C. L. Clum. Stochastic network analysis using moments. Unpublished M.S. thesis, The Ohio State Univ. (1983). 15. C. E. Clark. The greatest of a finite set of random variables. Ops Res. 0O, 145-162 (1961). 16. W. E. Wilhelm. Lognormal models of transient operations in the flexible manufacturing environment. J. Mtg Systems, forthcoming.
S, models of the maximum of several dependent random variables
345
17. N. L. Johnson and S. Kotz. Distributions in Statistics, Continuous Multivariate Distributions. Wiley, New York (1972). 18. J. Draper. Properties of a distribution resulting from certain simple transformations of the normal distribution. Biometrika 39, 290-301 (1952). 19. D. C. M. Leslie. Determination of parameters in the Johnson system of probability distributions. Biometrika 46,229-231 (1959). 20. J. R. Wilson. Fitting Johnson curves by univariate and multivariate data. Presented at the Winter Simulation Conf. (1983). 21. J. R. Wilson. Modeling multiattribute populations with Johnson's translation system. Technical Rep. Mech. Engng Dep., Univ. of Texas, Austin (1983).
APPENDICES
Derivations of the equations presented in Section 2 are outlined in the Appendices. Numbers in the left margin correspond to equation numbers in the text. These derivations parallel those used by Clark [15] to compute the moments of the maximum of bivariate normal random variables.
A p p e n d i x A : Pae (1) Product Moment E[YI Y2] = ~!~ ~l=sinh(n0 sinh(n2) d~N1.N2(nl,nz)dnldn 2 in which sinh(n 0 = (1/2)[exp(Ih +-O'XZl) e x p ( - ~ t l - o ~ z l ) ] and ~(.) is the density function of the standard normal. The second integral can be evaluated for K = + 1 as
I~(K) =
=!~exp[K(p,2 + 02z2)](1/(Zn)1~27/2) e x p [ ( - 1/2)(z2- p12zl) z]dz2
in which r = 1 - p t f . Upon completing the square I21K) = exp(K~t2 + r022/2) exp[KQ1202zd. Four terms must be integrated (for K = +1 and L = +1); each are of form
I( KL ) = exp( L~h + K~t2+ rozZ/2)(2~x)-1/2 ~l~exp{ (--1/2)[z lZ- 2z l( Kp1202 + LoO]dza and reduce to
I(KL) = e x p [ ( 1 / 2 ) ( ~ + ~ ) ] exp[K~tz+ L~h]exp[KLO,zOtO2] in which, it is noted, that only the last term is a function of P~z. Combining these 4 terms,
EI Y1Y2] = (1/4)exp[ (1/2)(o~2+ o~2)][S2exp(olo2o12)S3expolo2o12)]. Using this result and the identity Q(Y,,Y2) = (E[YIY2] -- E[Y~]E[Y,]E[Y2]}/V~/2(y,)v~/2(Y2), eqn. (1) may be derived directly by back-solving for ~12.
A p p e n d i x B: U i (6) Id¢ikdefined in eqn. (5) are derived using U
Wik =
i [sinh(~tk ~- OkZk)]i~)Zk'(Zk)
J d~Zk'lZk(Zk')dZk'dZk
(B1)
in which k = (1,2); k' = (1,2), k' v~ k
U = ~k ~- OkZk If sinh(nk) in eqn. (15) is expanded in its exponential form: [sinh(~tk
+ OkZk)] i ~
(l/2)i[exp(~tk + OkZk)-'exp(--~k
-- OkZk)] i
the terms defining Wik result. (7) Equation (15) differs from eqn. (4) in Wilhelm [161 only by the transformation applied to nk; results of the later case can be applied with i > IJ to obtain Mik in eqn. (3) or with i < 0 to obtain Mik. (2) Equation (2) results from using the integrand [~k -~- kk sinh(nk)] i
(B2)
U = ~tk + OkZk + ln(Xk/Xk').
(B3)
in eqn. (15) and upper limit of integration Equation (17) represents a special case of the general form U = sinh l[{~.ksinh(~t k + OkZk) -~ ~kt}/~.k ']
(B4)
which should be used. This general form is intractable, but may be evaluated for the special case in which ~l ~_2 and (~tk+OkZk) > 1 which, according to the following analysis, yields eqn (17). For k = 1 if (~t~ + otz~) > 1 then sinh(~h + o~z~) = (1/2) exp(~h + olzl) and if~t = ~_2 then U = sinh-l[(kJ2k2) exp(txl + olzl)]
346
W . E . Wtt.~tEt,m and if the bracketed term [(~.1/2~.z) exp(I,h +
ojzl)] is > 4
then U = en[(~.t/2k2) exp(p,1 + OIZl)], since sinh-I(Y) = fn[Y + (YZ+l)l/z] ~ fn[2Y] if Y > 4. This restriction requires IX, > In (8kJ~.~) + Okzk.
Appendix C: E[X 3 •max(X1, X2)] (9) Equation (9) may be derived with U = IXk + okz~, in U
Ek = f sinh(nk) OZk(Zk) l ~Zk'LZk(3k ') f sinh(n3)°z31zkzk,(z3)dZ3dZk'dZk
(Cl)
Each sinh(.) term gives rise to 2 exponential terms; their product, to 4 terms. Eight terms result ultimately since k must be evaluated at both 1 and 2. Each of these 8 terms is integrated using appropriate expressions needed to complete the square; eqn. (9) results. (8) Equation (8) may he derived using {(~.k/2)[exp(IXk + OkZk)--exp(--IXk -- OkZk)] + ~k} i
(C2)
with i = 1 and k = 1 or 2 as the first integrand in eqn. (19) and with k = 3 in the third integrand of eqn. (19). U should be defined according to eqn. (17).
0360-8352/86 $3.00 + 0.00 Pergamon Journals Ltd
Printed in Great Britain
Su M O D E L S OF T H E M A X I M U M OF S E V E R A L D E P E N D E N T RANDOM VARIABLES: DEVELOPMENT AND COMPUTER APPLICATION W . E.WILHELM Department of Industrial and Systems Engineering, The Ohio State University, Columbus OH 43210, U.S.A.
(Received 12 June 1985; received for publication 10 October 1985) Abstract - - An approach for modeling the distribution of the random variable (U) that is the
maximum of several dependent random variables is described in this paper. Applications of the model to various problems associated with the Industrial Engineering profession are described and one - - the accumulation of components in small-lot assembly systems - - is studied in detail. Numerical tests investigate the accuracy of the approach and identify several fundamental characteristics of the accumulation process. Apparent advantages of the model are that it requires computer run time that is a fraction of that required by Monte Carlo sampling (simulation) procedures and that it is readily amenable to parallel processing.
1. INTRODUCTION
A number of models in the field of Industrial Engineering require definition of the probability density function of the random variable that is the maximum of several dependent random variables. PERT networks, for example, require such a maximum to define the "realization" time of a node, the time at which all immediately preceding tasks are completed. According to Whitehouse[1], GERT network analysis entails the maximum operator to model AND nodes. More recently, Wilhelm and Ahmadi-Marandi[2] described an approach, based on the maximum operator, for modeling the transient performance of production/assembly queuing networks. Operation finish times are assumed to be related by the multivariate normal distribution and operation start times are then defined by evaluating the maximum of several dependent, normal random variables: Operation Start Time = max[X1, X2, X3] in which X l = finish time at the upstream station that readies the part, X2 = finish time at the upstream station that readies the subassembly, X 3 -- finish time of the operation that readies the assembly station. A similar evaluation has been used to estimate lifetimes of reliability systems in which components exhibit normally distributed lifetimes: System Lifetime = min [Xl, X2] = - m a x [ - X l , -X2] in which Xl and X2 are lifetimes of components 1 and 2 (resp.). Motivated by these applications, the objective of this study was to investigate the computational advantages of the S, distribution as a model of the random variable that is the maximum of several dependent random variables. The Su distribution, developed by Johnson[3], was selected for 2 reasons. First, the distribution is related to both the normal and lognormal distributions, since it is based on a transformation of a normally distributed random variable. Secondly, the distribution can be fit to a broad range of cases; in fact, it can fit all ([31, [52) combinations (i.e. measures of skewness and kurtosis) which lie in the half-plane below the "line" of ([51, [32) values that can be fit by the ~A~ ~0:4-E
335
336
W.E. WIHIEI,M
lognormal distribution. Based on the theory of extreme value distributions (e.g. see Johnson and Kotz[4]) it is expected that this region is the one of primary interest in practical applications. Prior theoretical research by mathematical statisticians (which is summarized by David[5], Galambos[6] and Gumbel[7]) has led to description of families of extreme value distributions, the limiting distributions (as J---> 2) of the maximum of J i.i.d. random variables. Most of the related, applied research has been motivated by the development of PERT and is summarized by Elmaghraby[8]. A variety of modeling approaches have been proposed including polynomial approximations by Martin[9], bounding distributions by Kleindorfer[10] and Shogan[ll], and hybrid schemes combining numerical and sampling techniques by Hartley and Wortham[12] and Ringer[13]. Clum[14] has recently investigated the use of Burr and S-D distributions to analyse PERT networks. None of these approaches offers the combination of accuracy and computational speed required in the applications described above. Clark[15] developed an alternate approach using the normal distribution as a modeling convention. He derived the (exact) first 4 moments of the random variable which is the maximum of 2 random variables that are related by the bivariate normal distribution. His convention is based on the assumption that this maximum is normally distributed so that the maximum of J dependent random variables can be modeled using a ( J - 1) stage (i.e. pairwise calculation indicated by max[Xl, X2 . . . . . XI] = m a x [ X j _ l . . . . max(X3, max(X2, Xl) • • .)]. This convention was applied to model assembly systems by Wilhelm and Ahmadi-Marandi[2] and was found to give "reasonable" estimates of system performance. Nevertheless, it is expected that accuracy can be improved by using conventions that permit more faithful representation of the maximum by a positively skewed distribution that more closely approximates third and fourth moments than does the (symmetrical) normal distribution. Wilhelm[16] developed and tested a convention based on the lognormai distribution, which is, of course, skewed and was shown to give accurate results with computation time a fraction of that required by sampling procedures. Results of the current research extend this line of investigation to define models based on the S, distribution. It is expected that this new convention will provide a fitting capability important in many applications. The primary advantage of these conventions is that they facilitate computer modeling of events that are defined by the maximum operator. Simulation can be applied to model these events, but runtime may be excessive and results may require involved statistical interpretation. Run times to implement the conventions are typically a fraction of the times that equivalent simulation models require; this attribute is particularly important for effective use of personal computers. In addition, the conventions can be implemented using parallel processors much more readily than simulation approaches. Finally, conventions offer convenient conceptual models which can lead to analytical models to replace simulation as the only available means of modeling certain processes. The body of this paper is organized in 3 sections. The S, convention, presented in Section 2, may be used to approximate the moments of the maximum of multivariate S, random variables, max[Xi, X2 . . . . . Xs], by taking the variables pairwise: U l = X 1
Ui = max[Xj, Uj-l] forj = 2, 3 . . . . .
J.
Equations which state the moments of Uj are given in Section 2 along with the product moment E[X3, U2] which is used to calculate the correlation r[Xj, Uj-*] required to estimate the moments of U. In addition, a procedure for fitting parameters to U~ is given. The convention is demonstrated in Section 3 in an application which involves coordinating the accumulation of components for assembly. This is an important practical problem which determines the amount of in-process inventory as well as schedule (due
S. models of the maximum of several dependent random variables
337
date) performance that can be achieved by an assembly system. Conclusions are given in Section 4. It should be noted that all derivations are outlined in the Appendices. 2. THE S, CONVENTION The which Z ~,, 6 Y
S, distribution is defined, according to Johnson [3], by Z = 5' + 6 s i n h - l ( Y ) , in = a unit normal variate = parameters of the S, distribution = the random variable described by the S, distribution
Equivalently, Z may be translated (i.e. N = ~t + o Z ) to yield the normally distributed random variable N with mean ~t and variance o 2, which, in turn, may be transformed to obtain Y, Y = sinh (N) = (e N - e-N)/2; the corresponding inverse transformation is N = sinh-l(Y) = f n [ Y + (y2 + 1)1/2]. The random variable Y that is defined by the transformation is said to have the 2-parameter S, distribution. The associated normal distribution, related by the transformation sinh(N), has parameters ~t = -5'/6 and o 2 = 1/6 z. If YI and Y2 have the bivariate S, distribution with coefficient of correlation Q(Y1, Y2), the correlation between the associated normal random variables Pl2 can be shown (see Appendix A) to be: Q,z = (1/olo2) tn{[S1 + (Sl 2 + 4S2S3)1/2]/(2S2)}
(1)
in which V[-] and E[.] are used as the variance and expected value operators (resp.) and S, = {o(Y,,
Y2)v~/z[Y,IVm[Y2I + E[Y, IE[Y2I)4
e x p [ - ( o l 2 + 022)/2]
Sz = exp(tx2 + Ix1) + exp(-Ix2 - gl) $3 = e x p ( - g 2 + gl) + exp(g2 - ~q). Flexibility may be enhanced by using the 4-parameter S, distribution to allow fitting the first 4 moments of a distribution, Based on a translation of Y, the random variable X = kY + ~ has mean E[X] = kE[Y] + ~ and central moments vi(X) = kivi(Y) for i >/2. It is easy to show that X and Y h a v e the same ([31, [32) values and that o(Yl, Y2) = 9(X1, X2). The 2-parameter case results from setting ~l = ~ 2 = 0 and kt = ~.2 = 1. A convention based on the S, distribution offers an appealing flexibility, since parameters of this distribution can be set to cover all points in the ([31, [3z) plane that lie "below" the lognormal (L) "line" (Johnson and Kotz [17]). For y < 0 the S, distribution is skewed to the right and may, therefore, provide a good fit for the distribution of max(Xl, )(2, . • . , X:) in which the Xj{j I J = 1, 2 . . . . . J} have the multivariate S,, distribution.
Moments of the maximum of bivariate Su random variables The i tn moment about zero for U = max(X1, X2), Ui, may be calculated for the special case (see Appendix B) in which El = ~ 2 and ~tk > ln(8kk'/~.k) + ok Z~ (for k = 1 and 2). Zk in the latter restriction may be replaced by 5, since Pr[ I Zk I <5] is close to 1.0 (for k = 1 and 2). These restrictions tend to limit the inherent flexibility of the 4-parameter distribution but may be satisfied in a number of situations of practical interest. In applications for which the restriction is reasonable, {/)i ] i = 1, 2, 3, 4} are defined (see Appendix B) by 2
0:
}-" [k,I/Clk + ~,F,] k:l 2
02 = 2 k=l
[kkW2k + 2~kkkl'9~k + ~2fk]
338
W.E. WILHELM 2
[kkW3k + 3~k;~kW2k + 3Xk ~kWlk
. =
+
~3fk ]
k--I 2
[)"kW4k + 4~k~'kW~k + 6~k~'kW2e +
=
kfk]" (2)
k=l
Terms in eqn (2) are defined in the following paragraphs. First, f,
=
-
.z
+
en(XdX2)l/a}
in which = distribution function of the standard normal a = [012 "Jr- 022 -
2~12G10"2] 1/2.
Due to symmetry, f2 =
in which the operator Q[.] is used to indicate that f2 can be obtained from fl by changing all subscripts 1(2) to 2(1) and will be used to denote a number of symmetrical relationships throughout this paper. If, in fact, ~ = ~2, 2 k=l
and the expressions for/..-/i in eqn (2) can be simplified appropriately. In order to define the l,~'ik terms in eqn (2), it is necessary to compute 4 additional terms, 3;/n, a;/il, a;/i2 and a5/i2: /~'/il = /-7'1 (I){[~1 -- ~'2 q-
iol (ol
-
Q12(J2) q'-
fn(~q/~.2)]/a},
(3)
in which /7,l = exp(ibtl + i2o12/2 ),
(4)
a;/il is/1;/n with i replaced by - i in both eqns (3) and (4), a;/i2 = f2[Ma], a;/i2 = fa[a;/i,] and L2 = f2[£.]. Using these terms, l'~'i~ (i.e. 14'ik for k = 1) is
l,~'i, = 2-i[1(4i, + (-1)iMi, + t,,] for i = 1, 2, 3, 4
(5)
in which
tll = O, t21 = - 2 , tBz = -614/11 and tal = -161'f'21 - 2; and, due to symmetry, 14'~2 = f~[W,I and ti2 = Q[t~a] for i = l, 2, 3, 4. It is interesting to note several relationships implicit in eqns (2)-(5). First, in the 2 parameter case [i.e. with 21 = 22 = 1.0 in eqn (3)] the ith m o m e n t about zero, if'i, of W = max((Y1, Y2) is
l'f'i = l'Vi, + I~i2.
(6)
Secondly, the ith moments about zero, M i, of M = max(L1, L2), in which L~ and L 2 a r e bivariate lognormals, is (see Wilhelm[16]) m i ~-- J~/il -{'- J~/i2
(7)
[still using ~-I = ~,2 = 1.0 in eqn (3)]. Finally, £1 (£2) in eqn (4) defines the i th moments about zero for £ | (£2), assuming 0l = 02 = 0.0 (see Johnson and Kotz[4].
S, models of the maximumof several dependent random variables
339
To summarize, the moments of U = max(X1, )(2) may be determined by computing, in order, (i) (ii) (iii) (iv) (v)
P12 using eqn (1), £1 and/]2 using eqn (4), Mil, ~/i2, /~/il a n d / ~ / i 2 using eqn (3), I~il and ff'n using eqn (5), Ui (for i = 1, 2, 3, 4) using eqn (2).
Product moment Specification of the Su convention is completed by defining the product moment E[X3-max(Xl, X2)] = E[X3. U] (see Appendix C): E[X3.max(X1, X2)] = Tit + Ti2
(for i = 1 only)
in which Til = ~-I~3V[/il jr_ ~l~3fl
+ (~.,K3)(D/4)[E,(1,1) + E l ( l , - 1 ) + E l ( - 1 , 1 ) + E , ( - 1 , - 1 ) ] + (~,k3)(O/2)[/~l(1,1) + / ~ ( 1 , - 1 ) + /~,(-1,1) + E , ( - 1 , - 1 ) ]
(8)
Ti2 = •[Til] and terms are defined in the following paragraphs. First, certain functions of the correlations are defined as R = 1 - P122 D = [(1 - ~122 - QI32 - Q232 q- 2QlaOI3Q23)/2R] G1 = (Q13 - ~12Q23)/R B12 = ol -- Q)1202 and G2 = ffaiG,], B2, = ffa[B,2]. The basic function Ek(',') is:
E , ( K , L ) = K L exp[IK~t3 + Lbt, + (P,a/2a 2) + QI] @[(~q - ~t2 + P,)/al in which, f o r K =
+landL
= _+1,
I=1 Al = IL01 + Ko3(Gt + (112G2) P1 = A1BI2 - GeKR0203 + ~n(~.i/~.2) Q1 = exp{(R/2a2)[Al2022 + BI22G22032 + 2A1Bt2G2K°2°3] en(?,,/?,2)[2P, - ln(~.,/~.2)]/R} -
By symmetry, A2 = Q[A,], P2 = ~2[Pl], Q2 = ~[Q,] and E2(K,L) = ffa[E,(K,L)]. It is of interest to note that Ee defines the product moment in the 2 parameter S,, case (with ~.l = ~.2 = 1.0): E[Ys.max(Yl, I12)] = E[Y3.W] = (D/4)[E,(1,1) + E , ( 1 , - 1 ) + E , ( - 1 , 1 ) + E , ( - 1 , - 1 ) + E2(1,1) + E2(1,-1) + E2(-1,1) + E 2 ( - 1 , - 1 ) ] .
(9)
Given these definitions,
E1(K,L) is E I ( K , L ) with I = 0 a n d / ) 2 ( K , L ) = f2[/~,(K,L)]. Parameter fitting procedure Implementation of the Su convention entails calculating the moments of U and then fitting the parameters of an S,, distribution to these moments. S, fitting procedures have
340
w.E. WII.HECM
been developed by Draper [18], Leslie [19], Johnson and Kotz [4], and more recently by Wilson [20, 21], but none of these is useful in the current application. A new procedure is, therefore, described in this subsection. If ~ and k are fixed (e.g. as they are in the 2 parameter case), the parameters (~t, o) of U may be determined by equating the moments of U [/_71 and U2m defined in eqn (2)] to those of an S. distribution: E[X ~] = /)~ in which,
E[X] = (k/2)e °-~/2[e ~' - e - " l + ~
(10)
E[X 2] = 0~2/4)[e 4°2/2 (e 2~t -- e -2~t) -- 2] + ~2 + ~.eo2/Z[e~t _ e - ' ] .
(11)
and
Equations (10) and (11) may be solved simultaneously for Su parameters ~t = In(F) and o 2 = 2 ln(r)
(12)
in which
~"
=
F 1/2, F = [~'I(F-1)][2(/21 - ~)/k] 2
and F is determined as a root of the polynomial equation F4-vl F3+6F 2-rlF+
1 =0
(13)
in which Vl = 4 + [2(/)1 - ~)/k]4[(4/)2/k 2) - (4~/L2)(2/)1 - ~) +
2]-'.
In the 4 parameters case ~ is fixed by the requirement that ~ = El = ~2, but parameters (k, tl, o) may be fitted to the first 3 moments of U, {/)i [ i = 1, 2, 3} defined by eqn (2). The third moment of X is E[X 3] = (~.3/4)e9°~-/2sinh(3g) + (~2 _ k2/4)(3?~eO-V2)sinh(g)
+ (3~k2/2)e 4''-/2 cosh(2~t) + (~2 _ 3~2/2)~,
(14)
so that (k, tx, o) may be estimated by solving eqns (10), (11), and (14) simultaneously. This approach requires solution of a high order polynomial and is, therefore, replaced by the algorithm stated below. (vi) S e t ~ = 1, k t = k m , B ~ = ~tm, Oe = o m , e > O . (vii) D ~= dE[X3]/dk = (0.75)(~f)2exp[4.5(ot)2]sinh[3~t ~] + (1.5)(k~)[3~ - (2 + /)1)]. = (viii) k t + l = k ~ + Ak. (ix) (x) Use ~ and U +1 in eqn. (12) to calculate ~V+1 and C +~. (xi) Use ~, U +1, I,d +1 and C +l in eqn (14) to calculate E[X3] t+l. (xii) If (1 - e)/)3 -< E[X3] c+1 < 1)3 go to step (viii). Else, £ <-- ¢J + 1 and go to step (vii). (xiii) Fitted parameters are e, )~+l, f + t , C + t . Calculate {E[Xi] I i = 1, 2, 3} and stop. It is implicitly assumed that the parameters of U are not dramatically different from those of Xm in which m is the index (i.e. 1 or 2) of max(~tl, ~t2). To begin, the parameters of U are set equal to those of variable m. Equation (14) is then used to improve the current estimate of ~., U. The rate of change of E[X 3] as a function of ~. (using current values ~tt, or), D t, is computed at step (vii) and used in step (viii) to estimate, according to a first-order Taylor series expansion, the increment Ak that would set E [ X 3] = 03 as desired. ~t' is updated using this increment [step (ix)] and then used, along with ~, to update estimates of (~t, o) using eqn (12). Revised parameters are used in eqn (14) to compute the corresponding value of/~[X3]. If this value is within E of U3, the procedure
S. models of the maximum of several dependent random variables
341
terminates; otherwise, it will iteratively improve parameter estimates until E[X1] = 0~, E[X2] = 02, and E[X3] is acceptably close to/-)3. 3. EXAMPLE APPLICATION
Components required to assemble a product must be accumulated before assembly can begin. For example, a printer circuit board cannot be assembled until the board and all electronic components have been accumulated; i.e. "kitted". Some products such as automobiles may require a kit for each subassembly as well as an accumulation of subassemblies before initiating final assembly. It is well known that the bulk of manufacturing is currently done in small lots to satisfy individual customer preferences. In small-lot assembly systems, components may be specially ordered rather than inventoried, because requirements are unique and products themselves are special ordered. Coordinating the input of materials from vendors and from production facilities is thus an important material management problem in small-lot assembly systems. By starting with the due date of a lot of assembled products and "backing off" an appropriate duration, the time at which an accumulation will be needed to initiate assembly can be readily determined. Material managers must then set a due date for the delivery of each component from its supplier. Due dates are thus used to coordinate material flows from vendors and production facilities to assembly operations. Since delivery time is a random variable in actual systems, component due dates must be set to hedge against the dual risks of early and late delivery. If a component is delivered well ahead of the time at which assembly is to begin, schedule performance will be good, but in-process inventory will be high. Conversely, if a component is delivered late, in-process inventory will be low, but schedule performance will be poor. Due dates must, therefore, be set to provide an appropriate balance between the risks of early and late delivery. The S, convention may be applied to study the accumulation process by using the model, with U ~ = X~ and Uj = max[U j - l , Xj] f o r j = 2, 3 . . . .
,J
in which X~ = delivery time of component j (j = 1, 2 . . . . . J) Uj = time at which component accumulation is complete. Computation procedures given previously in this paper may be used in the following manner to estimate the moments of accumulation time: (xiv) Set j = 1; (xv) Use procedure (i)-(v) to calculate moments of/_fl; (xvi) Use eqn (8) and standard indentities to calculate correlations: r[Xk, Uj] = r[Xk, max(U j - l , Xj] for k = j + l , . . . , J; (xvii) Use procedure (vi)-(xiii) to fit parameters of the S, distribution to moments calculated in step (xv) above. (xviii) I f j = J stop, moments of Uj have been calculated. Otherwise, s e t j = j + 1 and go to step (xv). Systematic use of step (xvi) yields pairwise correlations that are subsquently required in step (xv) to calculate moments. Applying the model in the context of component accumulation, E[Xi] is interpreted as the due date for component j and V[Xj] as the variance of delivery time. This model can be used to gain fundamental insights into characteristics of the accumulation process and could be extended, for example, to develop optimal due dates (although such work is beyond the scope of this paper).
342
W.E.
WII.ItELM
Pairwise correlations among delivery times (i.e. r[Xj, X / , ] ) are zero if the components are produced by independent vendors. However, if they are produced in a flexible manufacturing system (or some other common production facility) pairwise correlations will be nonzero. In an actual application, correlations could be estimated from historical experience or by modeling (e.g. see Wilhelm and Ahmadi-Marandi [2]). In this paper, 3 cases are considered for demonstration p u r p o s e s - Q(X1, X2) = 0.0, 0.15 and 0.30. The marginal mean and variance of component delivery times were assumed to equal 100 (~j = 0.0, ~.j -- 1.0, 5'j = -53.0683, and 6j = 10.0254) in all test cases reported in this section. These values were selected purposely, since they define a "worst case" in which it is most difficult to estimate the moments of Uj. If one delivery time had a mean that is much larger than other means, that delivery time would dominate the definition of accumulation time and make it easier to develop good estimates. Practical experience indicates that it is more difficult to meet assembly schedules if the number of components required in the accumulation J is large. Characteristics of the accumulation process are therefore investigated for 5 levels of J: 2, 3, 4, 5 and 10. In all, 15 cases were tested (3 levels of correlation for each of the 5 levels of J) by both the 2 and (restricted) 4 parameter S , distributions. The S, convention is based on the assumption that the Xj are related by the multivariate Su distribution; it imposes the approximation that Ui are also multivariate S, distributed. This introduces a potential source of error, because the true distribution of the maximum is not known. Since the S,, distribution can, in principle, be fit to a broad range of moments, the error may be of acceptable magnitude. Accuracy of the approach was, therefore, studied by comparing test results with estimates developed using a simulation model. Analytical models of the accumulation process are not available. Comparisons with the simulation model were made using several measures: -
-
KS -- KS test statistic related to the null hypothesis that the distribution function
developed from the simulation was that of the S, distribution (critical value for 30 cells and a = 0.50 is 0.241). - - A = (simulated average of E[UJ]) - (recursion estimate of E [ U ] ) . t = t-test statistic related to H0: A = 0, HA: A 4: O. _ _ ~2 = ~2_test statistic related to Ho: (simulated V[UJ]) = (recursion estimate of V[UJ] HA: (simulated V[UJ]) 4: (recursion estimate of V[UJ]). - - Computer runtimes (all tests were made on a PDP 11/44). -
-
Table 1. S,, convention
0
Case J
Average
Simulation* SD CPU (sec)
~t = ~2 = 0.0
Mean
~,1 = ~,2 = 1,0
Approximation SD CPU (sec)
KS
Results A
0.00
2 3 4 5 10
105.582 108.443 110.626 112.059 116.166
8.61 7.86 7.48 7.10 6.43
142.22 249.95 272.78 445.62 422.92
11)5.623 1118.578 110.549 112.001 116.155
8.74 8.04 7.58 7.24 6.29
0.150 0.285 0.751 1.500 2.784
0.014 0.028 0.018 0.022 0.031
-0.041 -0.142 0.077 0.058 11.1)11
0.15
2 3 4 5 10
105.125 107.712 109.696 111.073 114.859
.8.85 8.31 7.90 7.67 7.31
250.73 215.60 302.68 364.96 665.50
105.183 107.901 109.714 111.063 115.010
8.99 8.46 8,11 7.86 7.11
0.651 0.299 0.582 1.481 4.868
0.020 0.023 0.023 0.036 0.031
-0.1158 -0.190 -11.018 0.010 -0.151
0.30
2 3 4 5 10
104.619 106.923 108.696 110.002 113.482
9.06 8.75 8.32 8.25 8.16
151.21 190.23 264.37 384.02 555.26
104.71)2 107.157 108.828 110.091 113.850
9.23 8.84 8.58 8.38 7.78
0.151 0.351 0.466 0.683 2.800
0.016 0.022 0.025 0.033 0.032
-0.083 -I).234 -11.132 -0.089 -0.368
* 1000 replications. A = (Simulation average) - (approximation mean). KS = KS-test statistics based on 30 cells, S, distribution.
S,, m o d e l s o f t h e m a x i m u m o f s e v e r a l d e p e n d e n t r a n d o m v a r i a b l e s
343
Table 2 . 4 parameter S, convention (~l = ~2) Case J
Average
0.00
2 3 4 5 l0
105.582 108.443 110.626 112.059 116.166
8.61 7.86 7.48 7.10 6.43
142.22 249.95 272.78 445.62 422.92
105.623 108.585 110.549 112.001 116.153
8.73 8.(}4 7.58 7.24 6.28
0.265 0.649 0.817 1.683 3.831
0.014 0.028 0.018 0.022 0.029
-0.041 -0.142 (}.077 0.058 0.013
0.15
2 3 4 5 10
105.125 107.712 109.696 111.073 114.859
8.85 8.31 7.90 7.67 7.31
250.73 215.60 302.68 364.96 665.50
105.183 107.902 109.714 111.063 115.008
8.99 8.46 8.11 7.86 7.11
0.217 0.649 0.548 0.833 6.648
0.020 0.023 0.023 0.036 0.031
-0.058 -I}.190 -(}.1}18 0.010 -0.149
(}.30
2 3 4 5 10
104.619 106.923 11}8.696 I10.(X}2 113.482
9.06 8.75 8.32 8.25 8.16
151.21 190.23 264.37 384.02 555.26
104.702 107.157 108.827 110.090 113.848
9.23 8.84 8.58 8.38 7.78
0.151 0.351 0.435 0.583 3.850
0.016 0.022 0.025 0.033 0.032
-0.083 -0.234 -(}.131 -0.088 -0.366
0
Simulation* SD CPU (see)
Mean
Approximation SD CPU (sec)
KS
Results A
* 1000 replications. A = (Simulation average) - (approximation mean).
KS = KS-test statistics based on 30 cells, S. distribution.
Results for the 2 and (restricted) 4 parameter models are given in Tables 1 and 2 respectively. In all cases, the fitting procedure (step xvii) converged in one iteration (with = 0.001). None of the KS or t-statistics are significant (~ = 0.05) indicating that accumulation time, UJ, is approximately St, distributed and that the estimate of mean accumulation time, E[UJ], is acceptably close to the estimate devloped by the simulation model. All A values are small in absolute value, again indicating the accuracy with which the mean of accumulation time is estimated. X2 statistics are significant at an c~ = 0.05 level only for the variances V[UJ] in the 3 cases in which Q(X1, X2) = 0.30 and J = 10. Variance of accumulation time is thus estimated accurately in all but the most extreme cases tested. In all cases, runtime for the approximation model was a small fraction of that required for the simulation model to complete 1000 replications. The approximation model required approximately the time required to complete 5-7 simulation replications. In addition to these statements concerning the accuracy of the approximating procedure, it is possible to interpret results in relationship to the accumulation process. First, it is noted that a deterministic estimate which does not incorporate the variance of delivery time V[XI], significantly underestimates the time at which an accumulation will be complete. The deterministic estimate would be 100 in all cases and the true mean time of accumulation E[U] is seen to increase as the number of components J increases. Conversely, the variance of accumulation time decreases with an increase in J. This indicates that accumulations will be completed consistently at later times as J increases. Safety lead time, the difference between scheduled assembly time and component due date must increase as J increases to maintain an acceptable balance between the risks of early and late delivery. Results for different levels of pairwise correlations between component delivery times indicates that mean accumulation time decreases (and the variance of accumulation time increases) with an increase in pairwise correlation. This result suggests that assembly schedules may be more easily attained if components are produced in a common facility such as a flexible manufacturing system. This relationship has, apparently, not been reported previously in the literature and may suggest further investigation to understand advantages of integrated production/assembly facilities. 4. C O N C L U S I O N S
An approach based on the S,, distribution is presented in this paper as a means of modeling (approximately) the random variable that is the maximum of several dependCAIE
10:4-F
344
w . E . WILttEI.M
ent r a n d o m v a r i a b l e s . A p p l i c a t i o n s in which this a p p r o a c h might be used include m o d e l i n g t r a n s i e n t q u e u i n g n e t w o r k s (e.g. p r o d u c t i o n / a s s e m b l y s y s t e m s ) , P E R T a n d G E R T n e t w o r k s , a n d s y s t e m reliability. A p p l i c a t i o n to m o d e l the a c c u m u l a t i o n of c o m p o n e n t s in s m a l l - l o t a s s e m b l y s y s t e m s is s t u d i e d in detail. T h e a c c u m u l a t i o n p r o c e s s is i m p o r t a n t to m a t e r i a l flow m a n a g e m e n t , since it d e t e r m i n e s i n - p r o c e s s i n v e n t o r y levels as well as a s s e m b l y s c h e d u l e p e r f o r m a n c e . N u m e r i c a l tests i n d i c a t e t h a t the f u n d a m e n t a l a s s u m p t i o n u n d e r l y i n g the a p p r o a c h - that the m a x i m u m of m u l t i v a r i a t e S , r a n d o m v a r i a b l e s has the S , d i s t r i b u t i o n - - gives a g o o d a p p r o x i m a t i o n to the true d i s t r i b u t i o n . F u r t h e r m o r e , the a p p r o x i m a t i o n p r o c e d u r e gave e s t i m a t e s o f the m e a n of the m a x i m u m of d e p e n d e n t r a n d o m v a r i a b l e s which c o m p a r e f a v o r a b l y with t h o s e d e r i v e d f r o m a s i m u l a t i o n m o d e l . T h e s e results indicate that the a p p r o a c h gives a v i a b l e m o d e l for use in the a p p l i c a t i o n s d e s c r i b e d a b o v e . Test results also l e n d insight into the c h a r a c t e r i s t i c s of the a c c u m u l a t i o n p r o c e s s , i n d i c a t i n g that d e t e r m i n i s t i c p r o c e d u r e s for setting d u e d a t e s are i n a d e q u a t e a n d defining i m p o r t a n t r e l a t i o n s h i p s b e t w e e n the n u m b e r o f c o m p o n e n t s to be k i t t e d a n d the t i m e at which a c c u m u l a t i o n can be c o m p l e t e d . R e s u l t s also i n d i c a t e that a s s e m b l y s c h e d u l e s can be a t t a i n e d m o r e c o n s i s t e n t l y if c o m p o n e n t s are p r o d u c e d in a a c o m m o n facility. This p a p e r suggests f u r t h e r r e s e a r c h to d e v e l o p a l t e r n a t e a p p r o x i m a t i o n s for U as d e f i n e d in e q n (18) to relax r e s t r i c t i o n s which define a c c e p t a b l e v a l u e s o f rt a n d o. A d d i t i o n a l l y , e x t e n d e d m o d e l s of the a c c u m u l a t i o n p r o c e s s m a y be d e v e l o p e d to d e t e r m i n e o p t i m a l d u e d a t e s to b e t t e r c o o r d i n a t e m a t e r i a l flow in small-lot a s s e m b l y systems. P e r h a p s the m o s t significant c o n c l u s i o n s r e l a t e to c o m p u t e r a p p l i c a t i o n s . A l t h o u g h m a t h e m a t i c a l l y i n v o l v e d , the a p p r o a c h gives g o o d e s t i m a t e s of p e r f o r m a n c e with a r e l a t i v e l y low c o m p u t e r r u n t i m e . It is e x p e c t e d that this a d v a n t a g e will p r o m o t e use of s o m e m o d e l s o n p e r s o n a l c o m p u t e r s a n d s o l u t i o n of large-scale p r o b l e m s on larger c o m p u t e r s . G i v e n p a r a m e t e r v a l u e s , a n u m b e r o f e q u a t i o n s can be s o l v e d s i m u l t a n e o u s l y [e.g. steps (xv) a n d (xvi)] so the a p p r o a c h is well s u i t e d for p a r a l l e l processing. Simulation is t h e o n l y p r a c t i c a l a l t e r n a t i v e m e t h o d in m a n y a p p l i c a t i o n s , but it m a y r e q u i r e excessive r u n t i m e a n d is n o t r e a d i l y a d a p t e d to p a r a l l e l processing. Acknowledgements - - This material is based on work supported by the National Science Foundation under
grant no. MEA-8213196. The author gratefully acknowledges the work of Mr Surendra Saboo whose careful programming effort and thoughtful analysis of implementation requirements led to the numerical results reported in this paper. Also, constructive comments of the referees led to improvements in an earlier version of this paper and are hereby acknowledged.
REFERENCES 1. G. E. Whitehouse. Systems Analysis and Design Using Network Techniques. Prentice-Hall, Englewood Cliffs, NJ (1973). 2. W. E. Wilhelm and S. Ahmadi-Marandi. A methodology to describe operating characteristics of assembly systems, liE Trans. 3(14), 204-213 (1982). 3. N. L. Johnson. Systems on frequency curves generated by methods of translation. Biometrika 36, 149-176 (1949). 4. N. L. Johnson and S. Kotz. Distribution in Statistics: Univariate distributions - - 1. Wiley, New York (1970). 5. H. A. David. Order Statistics. Wiley, New York (1981). 6. J. Galambos. The Asymptotic Theory o f Extreme Order Statistics. Wiley, New York (1978). 7. E. J. Gumbel. Statistics o f Extremes. Columbia Univ. Press, New York (1958). 8. S. E. Elmaghraby. Activity Networks. Wiley, New York (1977). 9. J. J. Martin. Distribution of time through a directed acyclic network. Ops Res. 13(1), 46-66 (1965). 10. G. B. Kleindorfer. Bounding distributions for stochastic acyclic networks. Ops Res. 19, 1586-1601 (1971). 11. A. W. Shogan. Bounding distributions for a stochastic PERT network. Networks 7,359-381 (1977). 12. H. O. Hartley and A. W. Wortham. A statistical theory for PERT critical path analysis. Mgmt. Sci. 12(10), B469-B481 (1966). 13. L. J. Ringer. A statistical theory for PERT in which completion times of activities are inter-dependent. Mgmt. Sci. 17(11), 717-723 (1971). 14. C. L. Clum. Stochastic network analysis using moments. Unpublished M.S. thesis, The Ohio State Univ. (1983). 15. C. E. Clark. The greatest of a finite set of random variables. Ops Res. 0O, 145-162 (1961). 16. W. E. Wilhelm. Lognormal models of transient operations in the flexible manufacturing environment. J. Mtg Systems, forthcoming.
S, models of the maximum of several dependent random variables
345
17. N. L. Johnson and S. Kotz. Distributions in Statistics, Continuous Multivariate Distributions. Wiley, New York (1972). 18. J. Draper. Properties of a distribution resulting from certain simple transformations of the normal distribution. Biometrika 39, 290-301 (1952). 19. D. C. M. Leslie. Determination of parameters in the Johnson system of probability distributions. Biometrika 46,229-231 (1959). 20. J. R. Wilson. Fitting Johnson curves by univariate and multivariate data. Presented at the Winter Simulation Conf. (1983). 21. J. R. Wilson. Modeling multiattribute populations with Johnson's translation system. Technical Rep. Mech. Engng Dep., Univ. of Texas, Austin (1983).
APPENDICES
Derivations of the equations presented in Section 2 are outlined in the Appendices. Numbers in the left margin correspond to equation numbers in the text. These derivations parallel those used by Clark [15] to compute the moments of the maximum of bivariate normal random variables.
A p p e n d i x A : Pae (1) Product Moment E[YI Y2] = ~!~ ~l=sinh(n0 sinh(n2) d~N1.N2(nl,nz)dnldn 2 in which sinh(n 0 = (1/2)[exp(Ih +-O'XZl) e x p ( - ~ t l - o ~ z l ) ] and ~(.) is the density function of the standard normal. The second integral can be evaluated for K = + 1 as
I~(K) =
=!~exp[K(p,2 + 02z2)](1/(Zn)1~27/2) e x p [ ( - 1/2)(z2- p12zl) z]dz2
in which r = 1 - p t f . Upon completing the square I21K) = exp(K~t2 + r022/2) exp[KQ1202zd. Four terms must be integrated (for K = +1 and L = +1); each are of form
I( KL ) = exp( L~h + K~t2+ rozZ/2)(2~x)-1/2 ~l~exp{ (--1/2)[z lZ- 2z l( Kp1202 + LoO]dza and reduce to
I(KL) = e x p [ ( 1 / 2 ) ( ~ + ~ ) ] exp[K~tz+ L~h]exp[KLO,zOtO2] in which, it is noted, that only the last term is a function of P~z. Combining these 4 terms,
EI Y1Y2] = (1/4)exp[ (1/2)(o~2+ o~2)][S2exp(olo2o12)S3expolo2o12)]. Using this result and the identity Q(Y,,Y2) = (E[YIY2] -- E[Y~]E[Y,]E[Y2]}/V~/2(y,)v~/2(Y2), eqn. (1) may be derived directly by back-solving for ~12.
A p p e n d i x B: U i (6) Id¢ikdefined in eqn. (5) are derived using U
Wik =
i [sinh(~tk ~- OkZk)]i~)Zk'(Zk)
J d~Zk'lZk(Zk')dZk'dZk
(B1)
in which k = (1,2); k' = (1,2), k' v~ k
U = ~k ~- OkZk If sinh(nk) in eqn. (15) is expanded in its exponential form: [sinh(~tk
+ OkZk)] i ~
(l/2)i[exp(~tk + OkZk)-'exp(--~k
-- OkZk)] i
the terms defining Wik result. (7) Equation (15) differs from eqn. (4) in Wilhelm [161 only by the transformation applied to nk; results of the later case can be applied with i > IJ to obtain Mik in eqn. (3) or with i < 0 to obtain Mik. (2) Equation (2) results from using the integrand [~k -~- kk sinh(nk)] i
(B2)
U = ~tk + OkZk + ln(Xk/Xk').
(B3)
in eqn. (15) and upper limit of integration Equation (17) represents a special case of the general form U = sinh l[{~.ksinh(~t k + OkZk) -~ ~kt}/~.k ']
(B4)
which should be used. This general form is intractable, but may be evaluated for the special case in which ~l ~_2 and (~tk+OkZk) > 1 which, according to the following analysis, yields eqn (17). For k = 1 if (~t~ + otz~) > 1 then sinh(~h + o~z~) = (1/2) exp(~h + olzl) and if~t = ~_2 then U = sinh-l[(kJ2k2) exp(txl + olzl)]
346
W . E . Wtt.~tEt,m and if the bracketed term [(~.1/2~.z) exp(I,h +
ojzl)] is > 4
then U = en[(~.t/2k2) exp(p,1 + OIZl)], since sinh-I(Y) = fn[Y + (YZ+l)l/z] ~ fn[2Y] if Y > 4. This restriction requires IX, > In (8kJ~.~) + Okzk.
Appendix C: E[X 3 •max(X1, X2)] (9) Equation (9) may be derived with U = IXk + okz~, in U
Ek = f sinh(nk) OZk(Zk) l ~Zk'LZk(3k ') f sinh(n3)°z31zkzk,(z3)dZ3dZk'dZk
(Cl)
Each sinh(.) term gives rise to 2 exponential terms; their product, to 4 terms. Eight terms result ultimately since k must be evaluated at both 1 and 2. Each of these 8 terms is integrated using appropriate expressions needed to complete the square; eqn. (9) results. (8) Equation (8) may he derived using {(~.k/2)[exp(IXk + OkZk)--exp(--IXk -- OkZk)] + ~k} i
(C2)
with i = 1 and k = 1 or 2 as the first integrand in eqn. (19) and with k = 3 in the third integrand of eqn. (19). U should be defined according to eqn. (17).