- Email: [email protected]
Flowgraph models: a Bayesian case study in construction engineering
Journal of Statistical Planning and Inference 129 (2005) 181 – 193 www.elsevier.com/locate/jspi
Flowgraph models: a Bayesian case study in construction engineering Aparna V. Huzurbazar∗ Department of Mathematics and Statistics, University of New Mexico, Albuquerque, NM 87131-1141,USA Available online 21 August 2004
Abstract Flowgraph models are useful in a wide variety of systems engineering and survival analysis problems. They are especially useful for analyzing time-to-event data and constructing corresponding Bayes predictive distributions. When a continuous time semi-Markov process defines transition times between a finite number of states and interest focuses on estimating densities, survival/reliability and hazard functions, or predictive distributions, flowgraph models provide a way of presenting the model and an associated method for data analysis. The method is illustrated using data from a construction engineering project. © 2004 Elsevier B.V. All rights reserved. Keywords: Hazard; Predictive distribution; Reliability; Semi-Markov process; Time-to-event data
1. Introduction Flowgraph models are useful for modeling time to event data that result from a stochastic process. A flowgraph models potential outcomes, probabilities of outcomes, and waiting times for those outcomes to occur. Flowgraphs model semi-Markov processes and provide a practical alternative methodology for data analysis. They are useful when a continuous time semi-Markov process defines the transition times between states and interest focuses on estimating the density, reliability or survival function, or hazard function of the process. Given a stochastic process with conditionally independent states, a flowgraph model allows for the use of most standard waiting time distributions to model the different states. It provides ∗ Tel.: +1-505-277-2404; fax: +1-505-277-5505.
E-mail address: [email protected] (A.V. Huzurbazar). 0378-3758/$ - see front matter © 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.jspi.2004.06.046
182
A.V. Huzurbazar / Journal of Statistical Planning and Inference 129 (2005) 181 – 193
a method for accessing the waiting time distribution for any partial or total waiting time. Flowgraph models operate on moment generating functions (MGFs) and use saddlepoint approximations (cf. Daniels, 1954) to convert the MGFs to waiting time probability density functions (pdfs), cumulative distribution functions (CDFs), reliability or survival functions, and hazard functions. As a data analytic method, flowgraphs have distinct advantages over other methods for semi-Markov processes (cf. Ouhbi and Limnios, 1999, 2003). Flowgraphs handle censoring and can be used in either a frequentist or a Bayesian framework. Block diagrams and signal flowgraphs are widely used to represent engineering systems, especially in circuit analysis. Basic flowgraph ideas were developed in engineering, but they never incorporated probabilities, waiting times, or data analysis. The literature on flowgraph methods is vast. Introductions to flowgraph methods are contained in most circuit analysis or control systems textbooks such as Dorf and Bishop (1995), Gajic and Lelic (1996), and Whitehouse (1973). Statistical flowgraph models are based on flowgraph ideas but unlike their antecedents, flowgraph models can also be used to model and analyze data from complex stochastic systems. For literature on statistical flowgraph models see Butler and Huzurbazar (1997) and Huzurbazar (1999, 2000). This article focuses on an application of flowgraph modeling to a real data project involving planning in construction engineering. Section 2 provides a brief introduction to flowgraph models for engineering systems and Section 2.1 shows how to solve a flowgraph model. Section 3 describes the construction engineering project. Section 4 presents conclusions.
2. A brief introduction to flowgraph models for engineering systems Fig. 1 shows a complex system consisting of outcomes in series and cascaded in parallel with feedback loops. The system is an assembly line for a manufacturing process for car stereos. State 0 represents an initial detection of a problem with a stereo. The problem is categorized into one of two types of severity. If the severity is of type I, the system is in state 1 for repair of the item. Eventually, the problem is fixed and the item moves to state 3 where it is specifically inspected to make sure that the type I problem is fixed. If the problem is
Fig. 1. Flowgraph model for manufacturing system.
A.V. Huzurbazar / Journal of Statistical Planning and Inference 129 (2005) 181 – 193
183
not fixed, the item is returned to state 1, otherwise, it passes inspection and moves to state 5. Similarly, if the severity is of type II, the system is in state 2 for repair. Eventually, the problem is fixed and the item moves to state 4 where it is specifically inspected to make sure that the type II problem is fixed. If the problem is not fixed, the item is returned to state 2, otherwise, it passes inspection and moves to state 5. In a flowgraph model, the states or nodes represent outcomes. The nodes are connected by directed line segments called branches. These branches are labeled with transmittances. These transmittances are labeled with the “transition probability × moment generating function (MGF) the of waiting time distribution in the previous state” which is a quantity called the branch transmittance. In the figure, probabilities and moment generating functions of the waiting time distributions are shown as branch transmittances. The waiting times on the branches can be any parametric distributions that admit moment generating functions. Hence, the model is quite general in that Markovian assumptions are not made. We use the branch transmittances of a flowgraph model to solve for the MGF of the distribution of the waiting time of interest. This MGF is then converted to give a density, reliability, or hazard function. Quantities of interest include predicting the distribution of the total time until the item is up to standards, 0 → 5; predicting the waiting time in repair, say, 0 → 3 or 0 → 4, or 0 → 3 or 4; or predicting the total number of times an item fails inspection from a given state such as 3, i.e., the total number of times the transition 3 → 1 is made. More generally, flowgraph models can be used to assess system reliability and quantities such as the time to total and partial failure of a system. Huzurbazar (2000) gives details on using flowgraph models for computing the total and partial system failure for cellular telephone networks. Flowgraphs are distinct from other graphical models used in statistics. Directed graphical models use nodes to identify random variables and edges to model causal relationships. Flowgraphs, on the other hand, use directed graphs to model the outcomes of random variables. In a flowgraph, nodes identify the actual physical states of a system, edges model allowable transitions, probabilities of outcomes, and waiting times until the occurrence of outcomes. Complicated flowgraph systems require the use of reduction methods. This involves an algebraic technique called Mason’s rule developed in graph theory by Mason (1953) for solving systems of linear equations. While Mason’s rule does not involve probabilities or MGFs, it can be applied to statistical flowgraphs. The result from applying Mason’s rule is the MGF of the waiting time distribution of interest. The next section describes how to use Mason’s rule to solve the flowgraph of Fig. 1 and to convert the resulting MGF into a density, reliability, or hazard function. 2.1. Solving a flowgraph model Solving a flowgraph using Mason’s rule entails identifying all the paths and loops of a flowgraph, computing path and loop transmittances, and then solving for the overall MGF. Mason’s rule requires computing the transmittance for every distinct path from the initial state to the end state and adjusting for the transmittances of various loops. In Fig. 1 there are two paths from input to output, state 0 to state 5, given by 0 → 1 → 3 → 5 and 0 → 2 → 4 → 5 . Path transmittances are simply products. Fig. 1 contains two first-order loops. A first-order loop is any closed path that returns to the initiating node of the loop
184
A.V. Huzurbazar / Journal of Statistical Planning and Inference 129 (2005) 181 – 193
Table 1 Summary of paths and loops for the manufacturing system flowgraph model of Fig. 1 Paths
Transmittance
Path 1: 0 → 1 → 3 → 5 Path 2: 0 → 2 → 4 → 5
p01 p35 M01 (s)M13 (s)M35 (s) p02 p45 M02 (s)M24 (s)M45 (s)
1st-Order loops 1→3→1 2→4→2
p31 M13 (s)M31 (s) p42 M24 (s)M42 (s)
2nd-Order Loops 1 → 3 → 1 with 2 → 4 → 2
p31 p42 M13 (s)M31 (s)M24 (s)M42 (s)
Loops not touching path Path 1: 2 → 4 → 2 Path 2: 1 → 3 → 1
p42 M24 (s)M42 (s) p31 M13 (s)M31 (s)
without passing through any node more than once. In Fig. 1, these are 1 → 3 → 1 and 2 → 4 → 2. The transmittance of a first-order loop is the product of the individual branch transmittances involved in its passage. Table 1 gives the transmittances for these loops. Higher order loops are defined as follows. A jth-order loop consists of j non-touching firstorder loops. The transmittance of a higher-order loop is the product of the transmittances of the first-order loops it contains. Fig. 1 contains one second-order loop which is made up of the two non-touching first-order loops 1 → 3 → 1 and 2 → 4 → 2. There are no loops of order three or more. Table 1 gives a summary of the paths, loops, and their transmittances. The general form of Mason’s rule gives the MGF of the distribution of the overall waiting time, T, as M(s) =
1+
j i j (−1) Lj (s)] , j j (−1) Lj (s)
i Pi (s)[1 +
(2.1)
where Pi (s) is the transmittance for the ith path, Lj (s) in the denominator is the sum of the transmittances over the jth-order loops, and Lij (s) is the sum of the transmittances over jth-order loops sharing no common nodes with the ith path, i.e., loops not touching the path. Details of Mason’s rule can be found in D’Azzo and Houpis (1981) or Dorf and Bishop (1995). The overall MGF is also called the equivalent transmittance. Expression (2.1) for the MGF of the distribution T, the overall waiting time from 0 to 5, in Fig. 1 is M(s)=
p01 p35 M01 (s)M13 (s)M35 (s)[1−L11 (s)]+p02 p45 M02 (s)M24 (s)M45 (s)[1−L21 (s)] , 1 − L1 (s) + L2 (s)
(2.2)
A.V. Huzurbazar / Journal of Statistical Planning and Inference 129 (2005) 181 – 193
185
Table 2 Summary of waiting time distributions for the manufacturing system flowgraph of Fig. 1 Transmittance
Distribution
MGF
0→1
Inverse Gaussian (1 , 2 )
M01 (s) = exp 1 2 − 21 (22 − 2s)
1→3
Gamma (1 , 1 )
M13 (s) =
3→1
Exponential (1 )
M31 (s) =
3→5
Exponential (2 )
M35 (s) =
0→2
Exponential (3 )
M02 (s) =
2→4
Inverse Gaussian(2 , 2 )
M24 (s) = exp 1 2 − 21 (22 − 2s)
4→2
Gamma (2 , 2 )
M42 (s) =
4→5
Exponential (4 )
M45 (s) =
1 1 1 −s 1 1 −s 2 2 −s 3 3 −s
2 2 2 −s 4 4 −s
where L11 (s) = p42 M24 (s)M42 (s), L21 (s) = p31 M13 (s)M31 (s), L1 (s) = p31 M13 (s)M31 (s) + p42 M24 (s)M42 (s), L2 (s) = p31 M13 (s)M31 (s)p42 M24 (s)M42 (s). Note that partial waiting times, say from state 2 to state 5, can also be evaluated with these methods by solving for the required MGF, in this case, M25 (s). While the flowgraph model gives us the overall waiting time MGF, this MGF is not particularly useful unless it is converted to a quantity of interest such as a density, CDF, reliability, or hazard function. We use the manufacturing system of Fig. 1 for illustration. We assume a variety of waiting time distributions given in Table 2. If we had data on the real system, we would use that to select candidate distributions. Parameter estimation is done by standard methods such as maximum likelihood or the entire problem can be solved from the viewpoint of prediction using Bayesian methods as in Section 3. For now, we assume the following parameter values p01 = 1 − p02 = .6, p31 = 1 − p35 = .3, p42 = 1 − p45 = .5, 1 = 25, 2 = 3, 1 = 5, 1 = .8, 1 = .8, 2 = .5, 3 = .3 1 = 2, 2 = 4, 2 = 7, 2 = 3, 4 = 4. We have deliberately mixed distributions and chosen parameter values such that the transition 0 → 1 → 3 → 5 has much longer waiting time (heavier tailed overall distribution) than the competing path 0 → 2 → 4 → 5.
186
A.V. Huzurbazar / Journal of Statistical Planning and Inference 129 (2005) 181 – 193
The MGF for the overall waiting time using these distributions in (2.2) is 1 1 1 2 − 21 (22 −2s) M(s) = 0.42e −s 1 2 2 1 2 − 21 (22 −2.0s) (2 − s)−1 × 2 1 − 0.5e 2 − s 1 1 1 2 − 21 (22 −2s) −1 + 0.203 e 4 1 − 0.3 1 (1 − s) 1 − s × (3 − s)−1 (4 − s)−1
2 1 2 1 1 2 − 21 (22 −2s) −1 − log 1− 0.3 1(1 − s) − 0.5e 2 − s 1 − s 1 2 2 2 1 2 + 0.15 1 e1 2 − 1 (2 −2s) (1 − s)−1 . 1 − s 2 − s One method for converting MGFs into a densities is to use a saddlepoint approximation. A saddlepoint approximation for the density of one random variable T, is . (2.3) f (t) = {2K
(ˆs )}−1/2 exp{K(ˆs ) − sˆ t}, where K(·) = log M(·), cumulant generating function (CGF) for T , dK , K (·) = ds d 2K K
(·) = and ds 2
K (ˆs ) = t. The last equation, K (ˆs ) = t, is called the saddlepoint equation and sˆ is called the saddlepoint. Using this approximation requires only that the MGF of T exist in an open neighborhood of zero. The literature on saddlepoint approximations in statistics is extensive beginning with Daniels (1954). For an introduction to practical saddlepoint approximations for use in statistics see Huzurbazar, S. (1999). Note that our interest with flowgraph models is to convert the resulting MGF into a density, CDF, etc., and the saddlepoint approximation is one way to do this. There are numerical inversion techniques for Laplace transforms that could also be used. For example, see Ramaswami (2000). We use the saddlepoint with flowgraphs because it is an analytical approximation that is fast, accurate, and well-suited in this case. While expressions for flowgraph MGFs may look complicated, flowgraph algebra lends itself to symbolic computation. The flowgraph structure can be computed using a symbolic algebra package such as MAPLE, which in turn can be used to generate C++ or FORTRAN code for the necessary functions, K(·), K (·), and K
(·), that are inputs to the saddlepoint approximation. The resulting density and hazard functions are given in Fig. 2. The bump in the density is reflective of the manner in which our finite mixture was chosen. This type of modeling allows us to realistically model the shape of the hazard to include portions that are increasing,
A.V. Huzurbazar / Journal of Statistical Planning and Inference 129 (2005) 181 – 193
187
0.08
density
0.06 0.04 0.02 0.00 0
10
20 30 time (in days)
40
50
0
10
20 30 time (in days)
40
50
hazard
0.15 0.10 0.05 0.0
Fig. 2. Density and hazard functions for manufacturing system.
decreasing, or flat (cf. Aalen and Gjessing, 2001). For a given system, there is a wide scope for fine tuning the underlying distributions to reflect the actual system behavior.
3. Construction engineering application A construction engineering project examined the impact of communication and planning on performance (cf. Garcia, 1997). Data on every aspect of operations was collected for twenty firms in New Mexico, USA. Monthly performance is one of the main quantities of interest. Monthly performance is loosely defined as whether the project completed the month on schedule and whether it will complete the next month on schedule. To evaluate completion times, current methods used (reluctantly) by these firms are Critical Path Method (CPM) and Project Evaluation and Review Technique (PERT). Both of these methods are commonly used in management and engineering to identify potential problems in a project plan but it is accepted that they have many shortcomings. See Pritsker et al. (1989) and AbouRizk and Wales (1997) for details. They are used at the beginning of a project and are not very adaptable to the project once the work has begun. They do not provide much assistance to the project manager in terms of timely updates to the project schedule. See Leung and Tam (1999) for details. CPM is also limited to analyzing deterministic systems with no feedback or failures in the project. Projects are always assumed to complete successfully. PERT is limited to using beta and normal distributions. Some current trends include the
188
A.V. Huzurbazar / Journal of Statistical Planning and Inference 129 (2005) 181 – 193
Fig. 3. Flowgraph model for construction engineering project.
use of genetic algorithms (cf. Li et al., 1999). Genetic algorithms and other evolutionary programming methods require large amounts of data for the algorithm’s learning process and a great deal of computation to implement. We use these data to illustrate a Bayesian predictive analysis from the flowgraph model. Fig. 3 shows a flowgraph model for the construction engineering project. State 0 is the planning state and state 1 is the design phase. From state 1, transition is possible to state 6, modifications, or into the monthly planning phase, state 2. The loop from 1 → 6 → 1 represents the process of design with modifications. After the initial planning, the monthly phase begins, represented by the boxed subsystem in Fig. 3. The time frame for the subsystem is one month (=23 working days). State 2 is the monthly planning state, state 3 is the start of the monthly work, state 4 represents delays in the monthly work which can result in modifications to the main project, indicated by a line from 4 → 6, or eventually lead to completion, state 5. Analysis from a flowgraph model gives an entire waiting time distribution as well as the system reliability and hazard functions for any total or partial waiting time. For decision making, this immediately provides the mean, median, or any percentile of the waiting time distribution for a variety of input distributions. The analysis from the current month along with prior information about potential delays gives a predictor distribution for the next month. This gives the project team a figure for how many days behind schedule they are likely to be for the next month allowing them to make necessary modifications such as reorganizing workload or hiring additional workers. Other methods generally only provide the mean and median for only a few distributions. Our data for the twenty firms is restricted to states 2–5 and consist of observations only on the monthly phase so our analysis assumes that from state 4, only transitions to state 5 are possible (p45 = 1). With appropriate data, it would be possible to analyze the entire overall schedule with the flowgraph model given in Fig. 3. Six of the firms experienced no delays. For these firms, we observe transitions from MONTHLY PLAN, to START, to COMPLETION 2 → 3 → 5. Of the 14 firms that experienced delays, 10 were able to recover from these and complete the month on schedule, 4 → 5. The remaining 4 firms were irrecoverably delayed,
A.V. Huzurbazar / Journal of Statistical Planning and Inference 129 (2005) 181 – 193
189
Histogram for 3 > 4
0.12
0.10
0.08
0.06
0.04 Gamma 0.02
0.00
2
6
10
14
18
22
26
30
34
38
Waiting Time Fig. 4. Histogram for data on the 3 → 4.
indicated by the transition 4 → 6. We have no subsequent data for them on the current month and therefore, in the present analysis, we will not analyze this transition further. We use histograms of the data to suggest parametric models to be used with flowgraphs. For example, Fig. 4 presents a histogram of data on the 3 → 4 transition. For illustration, an estimated gamma density is overlaid in the figure. After examining histograms for each state-to-state transition of the data, we make the following parametric assumptions about the branch transmittances: • • • •
Branch 2 → 3, constant 4 days, MGF M23 (s) = exp(4s); Branch 3 → 4, gamma(, ) with mean /, MGF M34 (s) = [/( − s)] , for s < ; Branch 4 → 5, constant 18 days, MGF M45 (s) = exp(18s); Branch 3 → 5, inverse Gaussian (1 , 2 ) with density f (t | 1 , 2 )= √
1 exp{−(21 /t+22 t)/2+1 2 } for t > 0, 1 > 0, 2 0 2t 3/2
and MGF M35 (s) = exp{1 2 −
21 (22 − 2s)},
s<
22 . 2
190
A.V. Huzurbazar / Journal of Statistical Planning and Inference 129 (2005) 181 – 193
Data on 2 → 3 are all within a small fraction of 4 days and data on the 4 → 5 transition are clustered around 18 days. We treat these as roughly constant times to demonstrate how flowgraphs handle constant variables. The 3 → 5 transition has only 6 observations. For illustration, we use an inverse Gaussian distribution to allow for a heavy tail. This is motivated largely by the problem at hand. The 3 → 5 transition means that work is completed without delays. In the construction industry, this is rarely the case. The heavy tail allows for longer waiting times when delays exist but are not recorded. Using (2.1) the MGF of the overall waiting time distribution from 2 → 5. There are two paths through the system, P1 (s) = p34 M23 (s)M34 (s)M45 (s) P2 (s) = p35 M23 (s)M35 (s) where p34 + p35 = 1, and we take p12 = 1 and p45 = 1 for the monthly phase analysis. The MGF of the total waiting time distribution for monthly time to completion is M25 (s) = e p34 4s
−s
e
18s
2 2 + e (1 − p34 ) exp 1 2 − 1 (2 − 2s) (3.1) 4s
for s < min(, 22 /2), 1 > 0, 2 > 0, > 0, > 0, 0 p34 leq1. The model can accommodate any transition time distribution, discrete or continuous. In particular, if 4 → 5 took on the values 18 and 19 with probabilities 23 and 13 , the term e18s in (3.1) would be replaced by (2/3)e18s + (1/3)e19s . Our goal is to compute the Bayes predictive density of the overall waiting time from planning to completion, 2 → 5, for a given month. We let T represent this random time. The Bayes predictive density of a future observable Z given data D about a system with parameter vector is f (z | D) =
f (z | )L( | D)() d ≡ E( | D) {f (z | )}, L( | D)() d
(3.2)
where z has density f (z | ), L( | D) is the likelihood function, and () is the prior. Our parameter vector is = (p34 , , , 1 , 2 ). To compute the Bayes Predictive density of (3.6), we assign a U nif (0, 1) prior for p34 and we assume diffuse U nif (0, 20) priors on all of the other parameters. The likelihood is separable in p34 , (, ), and (1 , 2 ). This results in a beta(15, 7) posterior distribution on p34 and the posterior on the remaining parameters is given by
(, , 1 , 2 | D) ∝
()
14 14
ti−1
e−
14
i=1 ti
i=1 20 20 2 2 1 × 61 exp − 1 − 2 t j + 6 1 2 2 tj 2 j =15
j =15
(3.3)
A.V. Huzurbazar / Journal of Statistical Planning and Inference 129 (2005) 181 – 193
191
0.04
density
0.03
0.02
0.01
0.00 0
10
20
30
40
50
t
Fig. 5. Bayes predictive density function for the completion time of the construction engineering project.
for 0 < , , 1 , 2 < 20. where the ti , . . . , t14 are the observed waiting times for the transition 3 → 4 and t15 , . . . , t20 are the observed waiting times for the transition 3 → 5. Note that the waiting times for transitions 2 → 3 and 4 → 5 are constant. Monte Carlo methods are used to perform the integration in (3.2) using a random sample of values from the posterior in (3.3) and from the beta(15, 7) posterior. Posterior sampling is done via slice sampling (cf. Neal, 2003), a Markov chain Monte Carlo method. Computationally, this involves using slice sampling to get a random sample 1 , . . . , m from the posterior. For each k value, we evaluate fˆ(z | k ) and L(k | D). fˆ(z | k ) is computed using the saddlepoint approximation (2.3) to convert the underlying flowgraph model MGF (3.1) into a density value. Given data D, the likelihood function for each k is evaluated based on the parametric assumptions on the branch waiting times. Figs. 5 and 6 show the Bayes predictive density and hazard function for this project. The predictive mean is 32.29 days and the standard deviation is 141.91 days so we do not expect to remain on schedule. The large standard deviation is reflective of the very heavy tailed distribution for the 3 → 5 transition and the fact that the once the project falls behind schedule it takes much longer to complete due to compounding problems. The hazard increases steadily until about 30 days and then drops off. This is reflective of the fact that project delays make it hard to complete in the 20 additional days shown in the figure. However, once the problems are resolved, the the hazard for the completion time will start to increase again in the furture. The 25th, 50th, 75th, 90th, 95th, and 99th percentiles are: 20.99, 27.28, 34.25, 44.57, 55.33 and 103.76 days. Any other desired quantity is also easily obtained from the predictive density. Alternatively, we could have used maximum likelihood to estimate the density of Z as f (z | ˆ ). The flowgraph provides the MGF which can be evaluated at ˆ , the MLE, and inverted to give an approximate density, fˆ(z | ˆ ), using the saddlepoint approximation. The
192
A.V. Huzurbazar / Journal of Statistical Planning and Inference 129 (2005) 181 – 193
0.10
hazard
0.08
0.06
0.04
0.02
0.00 0
10
20
30
40
50
t
Fig. 6. Bayes predictive hazard function for the completion time of the construction engineering project.
shape of the hazard is increasing and decreasing at varying rates. This is typical of hazards from flowgraph models since we are modeling the actual waiting times and not making assumptions about the hazards.
4. Conclusion Flowgraphs provide a data-analytic approach for modeling time-to-event data. Computational aspects begin with determining a flowgraph model for a given system. Estimation based on the flowgraph model can be performed in the Bayesian framework or via maximum likelihood. MGFs of total waiting times of interest are then computed from the flowgraph model. Inversion methods such as saddlepoint approximations can be used to convert this into a density, CDF, reliability, and hazard functions. Flowgraph models allow data analysis for semi-Markov processes using non-exponential waiting times. They also allow analysis of systems with feedforward and feedback loops in the presence of general waiting time distributions. Flowgraphs model semi-Markov processes that are applicable to a wide variety of situations. In the engineering setting, they can be used for systems analysis and design of complex systems.
Acknowledgements The author thanks Sandra Garcia, Federal Highway Administration, Department of Transportation, Boise, ID, for providing the construction engineering data and two anonymous referees for helpful comments.
A.V. Huzurbazar / Journal of Statistical Planning and Inference 129 (2005) 181 – 193
193
References Aalen, O., Gjessing, H., 2001. Understanding the shape of the hazard rate: a process point of view. Statist. Sci. 16, 1–13. AbouRizk, S., Wales, R., 1997. Combined discrete-event/continuous simulation for project planning. J. Construct. Eng. Manage. 123, 11–20. D’Azzo, J., Houpis, C., 1981. Linear Control System Analysis and Design: Conventional and Modern, McGrawHill, New York. Butler, R., Huzurbazar, A., 1997. Stochastic network models for survival analysis. J. Amer. Statist. Assoc. 92, 246–257. Daniels, H., 1954. Saddlepoint approximations in statistics. Ann. Math. Statist. 25, 631–650. Dorf, R., Bishop, R., 1995. Modern Control Systems, Addison-Wesley, Reading, MA. Gajic, Z., Lelic, M., 1996. Modern Control Systems Engineering, Prentice-Hall, New York. Garcia, S., 1997. Planning and Communication in Construction: Impacts on Performance. M.S. Thesis, University of New Mexico. Huzurbazar, A., 1999. Flowgraph models for generalized phase type distributions with non-exponential waiting times. Scand. J. Statist. 26, 145–157. Huzurbazar, A., 2000. Modeling and analysis of engineering systems data using flowgraph models. Technometrics 42, 300–306. Huzurbazar, S., 1999. Practical saddlepoint approximations. Amer. Statist. 53, 225–232. Leung, A., Tam, C., 1999. Models for assessing hoisting times of tower cranes. J. Construct. Eng. Manage. 125, 385–391. Li, H., Cao, J., Love, P., 1999. Using machine learning and GA to solve time-cost trade-off problems. J. Construct. Eng. Manage. 125, 347–353. Mason, S.J., 1953. Feedback theory—some properties of signal flow graphs. Proc. Inst. Radio Eng. 41, 1144–1156. Neal, R., 2003. Slice Sampling. Ann. Statist. 31, 705–767. Ouhbi, B., Limnios, N., 1999. Nonparametric estimation for semi-Markov processes based on its hazard rate functions. Statist. Infer. Stoch. Process. 2, 151–173. Ouhbi, B., Limnios, N., 2003. Nonparametric reliability estimation of semi-Markov processes. J. Statist. Plann. Infer. 109, 155–165. Pritsker, A., Sigal, C., Hammesfahr, R., 1989. SLAM II Network Models for Decision Support, Prentice-Hall, New Jersey. Ramaswami, V., 2000. Algorithmic analysis of stochastic models: the changing face of mathematics. Ramanujan Endowment Lecture at Anna University, Chennai, India. Whitehouse, G., 1973. Systems Analysis and Design Using Network Techniques. Prentice Hall, NJ.
Flowgraph models: a Bayesian case study in construction engineering Aparna V. Huzurbazar∗ Department of Mathematics and Statistics, University of New Mexico, Albuquerque, NM 87131-1141,USA Available online 21 August 2004
Abstract Flowgraph models are useful in a wide variety of systems engineering and survival analysis problems. They are especially useful for analyzing time-to-event data and constructing corresponding Bayes predictive distributions. When a continuous time semi-Markov process defines transition times between a finite number of states and interest focuses on estimating densities, survival/reliability and hazard functions, or predictive distributions, flowgraph models provide a way of presenting the model and an associated method for data analysis. The method is illustrated using data from a construction engineering project. © 2004 Elsevier B.V. All rights reserved. Keywords: Hazard; Predictive distribution; Reliability; Semi-Markov process; Time-to-event data
1. Introduction Flowgraph models are useful for modeling time to event data that result from a stochastic process. A flowgraph models potential outcomes, probabilities of outcomes, and waiting times for those outcomes to occur. Flowgraphs model semi-Markov processes and provide a practical alternative methodology for data analysis. They are useful when a continuous time semi-Markov process defines the transition times between states and interest focuses on estimating the density, reliability or survival function, or hazard function of the process. Given a stochastic process with conditionally independent states, a flowgraph model allows for the use of most standard waiting time distributions to model the different states. It provides ∗ Tel.: +1-505-277-2404; fax: +1-505-277-5505.
E-mail address: [email protected] (A.V. Huzurbazar). 0378-3758/$ - see front matter © 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.jspi.2004.06.046
182
A.V. Huzurbazar / Journal of Statistical Planning and Inference 129 (2005) 181 – 193
a method for accessing the waiting time distribution for any partial or total waiting time. Flowgraph models operate on moment generating functions (MGFs) and use saddlepoint approximations (cf. Daniels, 1954) to convert the MGFs to waiting time probability density functions (pdfs), cumulative distribution functions (CDFs), reliability or survival functions, and hazard functions. As a data analytic method, flowgraphs have distinct advantages over other methods for semi-Markov processes (cf. Ouhbi and Limnios, 1999, 2003). Flowgraphs handle censoring and can be used in either a frequentist or a Bayesian framework. Block diagrams and signal flowgraphs are widely used to represent engineering systems, especially in circuit analysis. Basic flowgraph ideas were developed in engineering, but they never incorporated probabilities, waiting times, or data analysis. The literature on flowgraph methods is vast. Introductions to flowgraph methods are contained in most circuit analysis or control systems textbooks such as Dorf and Bishop (1995), Gajic and Lelic (1996), and Whitehouse (1973). Statistical flowgraph models are based on flowgraph ideas but unlike their antecedents, flowgraph models can also be used to model and analyze data from complex stochastic systems. For literature on statistical flowgraph models see Butler and Huzurbazar (1997) and Huzurbazar (1999, 2000). This article focuses on an application of flowgraph modeling to a real data project involving planning in construction engineering. Section 2 provides a brief introduction to flowgraph models for engineering systems and Section 2.1 shows how to solve a flowgraph model. Section 3 describes the construction engineering project. Section 4 presents conclusions.
2. A brief introduction to flowgraph models for engineering systems Fig. 1 shows a complex system consisting of outcomes in series and cascaded in parallel with feedback loops. The system is an assembly line for a manufacturing process for car stereos. State 0 represents an initial detection of a problem with a stereo. The problem is categorized into one of two types of severity. If the severity is of type I, the system is in state 1 for repair of the item. Eventually, the problem is fixed and the item moves to state 3 where it is specifically inspected to make sure that the type I problem is fixed. If the problem is
Fig. 1. Flowgraph model for manufacturing system.
A.V. Huzurbazar / Journal of Statistical Planning and Inference 129 (2005) 181 – 193
183
not fixed, the item is returned to state 1, otherwise, it passes inspection and moves to state 5. Similarly, if the severity is of type II, the system is in state 2 for repair. Eventually, the problem is fixed and the item moves to state 4 where it is specifically inspected to make sure that the type II problem is fixed. If the problem is not fixed, the item is returned to state 2, otherwise, it passes inspection and moves to state 5. In a flowgraph model, the states or nodes represent outcomes. The nodes are connected by directed line segments called branches. These branches are labeled with transmittances. These transmittances are labeled with the “transition probability × moment generating function (MGF) the of waiting time distribution in the previous state” which is a quantity called the branch transmittance. In the figure, probabilities and moment generating functions of the waiting time distributions are shown as branch transmittances. The waiting times on the branches can be any parametric distributions that admit moment generating functions. Hence, the model is quite general in that Markovian assumptions are not made. We use the branch transmittances of a flowgraph model to solve for the MGF of the distribution of the waiting time of interest. This MGF is then converted to give a density, reliability, or hazard function. Quantities of interest include predicting the distribution of the total time until the item is up to standards, 0 → 5; predicting the waiting time in repair, say, 0 → 3 or 0 → 4, or 0 → 3 or 4; or predicting the total number of times an item fails inspection from a given state such as 3, i.e., the total number of times the transition 3 → 1 is made. More generally, flowgraph models can be used to assess system reliability and quantities such as the time to total and partial failure of a system. Huzurbazar (2000) gives details on using flowgraph models for computing the total and partial system failure for cellular telephone networks. Flowgraphs are distinct from other graphical models used in statistics. Directed graphical models use nodes to identify random variables and edges to model causal relationships. Flowgraphs, on the other hand, use directed graphs to model the outcomes of random variables. In a flowgraph, nodes identify the actual physical states of a system, edges model allowable transitions, probabilities of outcomes, and waiting times until the occurrence of outcomes. Complicated flowgraph systems require the use of reduction methods. This involves an algebraic technique called Mason’s rule developed in graph theory by Mason (1953) for solving systems of linear equations. While Mason’s rule does not involve probabilities or MGFs, it can be applied to statistical flowgraphs. The result from applying Mason’s rule is the MGF of the waiting time distribution of interest. The next section describes how to use Mason’s rule to solve the flowgraph of Fig. 1 and to convert the resulting MGF into a density, reliability, or hazard function. 2.1. Solving a flowgraph model Solving a flowgraph using Mason’s rule entails identifying all the paths and loops of a flowgraph, computing path and loop transmittances, and then solving for the overall MGF. Mason’s rule requires computing the transmittance for every distinct path from the initial state to the end state and adjusting for the transmittances of various loops. In Fig. 1 there are two paths from input to output, state 0 to state 5, given by 0 → 1 → 3 → 5 and 0 → 2 → 4 → 5 . Path transmittances are simply products. Fig. 1 contains two first-order loops. A first-order loop is any closed path that returns to the initiating node of the loop
184
A.V. Huzurbazar / Journal of Statistical Planning and Inference 129 (2005) 181 – 193
Table 1 Summary of paths and loops for the manufacturing system flowgraph model of Fig. 1 Paths
Transmittance
Path 1: 0 → 1 → 3 → 5 Path 2: 0 → 2 → 4 → 5
p01 p35 M01 (s)M13 (s)M35 (s) p02 p45 M02 (s)M24 (s)M45 (s)
1st-Order loops 1→3→1 2→4→2
p31 M13 (s)M31 (s) p42 M24 (s)M42 (s)
2nd-Order Loops 1 → 3 → 1 with 2 → 4 → 2
p31 p42 M13 (s)M31 (s)M24 (s)M42 (s)
Loops not touching path Path 1: 2 → 4 → 2 Path 2: 1 → 3 → 1
p42 M24 (s)M42 (s) p31 M13 (s)M31 (s)
without passing through any node more than once. In Fig. 1, these are 1 → 3 → 1 and 2 → 4 → 2. The transmittance of a first-order loop is the product of the individual branch transmittances involved in its passage. Table 1 gives the transmittances for these loops. Higher order loops are defined as follows. A jth-order loop consists of j non-touching firstorder loops. The transmittance of a higher-order loop is the product of the transmittances of the first-order loops it contains. Fig. 1 contains one second-order loop which is made up of the two non-touching first-order loops 1 → 3 → 1 and 2 → 4 → 2. There are no loops of order three or more. Table 1 gives a summary of the paths, loops, and their transmittances. The general form of Mason’s rule gives the MGF of the distribution of the overall waiting time, T, as M(s) =
1+
j i j (−1) Lj (s)] , j j (−1) Lj (s)
i Pi (s)[1 +
(2.1)
where Pi (s) is the transmittance for the ith path, Lj (s) in the denominator is the sum of the transmittances over the jth-order loops, and Lij (s) is the sum of the transmittances over jth-order loops sharing no common nodes with the ith path, i.e., loops not touching the path. Details of Mason’s rule can be found in D’Azzo and Houpis (1981) or Dorf and Bishop (1995). The overall MGF is also called the equivalent transmittance. Expression (2.1) for the MGF of the distribution T, the overall waiting time from 0 to 5, in Fig. 1 is M(s)=
p01 p35 M01 (s)M13 (s)M35 (s)[1−L11 (s)]+p02 p45 M02 (s)M24 (s)M45 (s)[1−L21 (s)] , 1 − L1 (s) + L2 (s)
(2.2)
A.V. Huzurbazar / Journal of Statistical Planning and Inference 129 (2005) 181 – 193
185
Table 2 Summary of waiting time distributions for the manufacturing system flowgraph of Fig. 1 Transmittance
Distribution
MGF
0→1
Inverse Gaussian (1 , 2 )
M01 (s) = exp 1 2 − 21 (22 − 2s)
1→3
Gamma (1 , 1 )
M13 (s) =
3→1
Exponential (1 )
M31 (s) =
3→5
Exponential (2 )
M35 (s) =
0→2
Exponential (3 )
M02 (s) =
2→4
Inverse Gaussian(2 , 2 )
M24 (s) = exp 1 2 − 21 (22 − 2s)
4→2
Gamma (2 , 2 )
M42 (s) =
4→5
Exponential (4 )
M45 (s) =
1 1 1 −s 1 1 −s 2 2 −s 3 3 −s
2 2 2 −s 4 4 −s
where L11 (s) = p42 M24 (s)M42 (s), L21 (s) = p31 M13 (s)M31 (s), L1 (s) = p31 M13 (s)M31 (s) + p42 M24 (s)M42 (s), L2 (s) = p31 M13 (s)M31 (s)p42 M24 (s)M42 (s). Note that partial waiting times, say from state 2 to state 5, can also be evaluated with these methods by solving for the required MGF, in this case, M25 (s). While the flowgraph model gives us the overall waiting time MGF, this MGF is not particularly useful unless it is converted to a quantity of interest such as a density, CDF, reliability, or hazard function. We use the manufacturing system of Fig. 1 for illustration. We assume a variety of waiting time distributions given in Table 2. If we had data on the real system, we would use that to select candidate distributions. Parameter estimation is done by standard methods such as maximum likelihood or the entire problem can be solved from the viewpoint of prediction using Bayesian methods as in Section 3. For now, we assume the following parameter values p01 = 1 − p02 = .6, p31 = 1 − p35 = .3, p42 = 1 − p45 = .5, 1 = 25, 2 = 3, 1 = 5, 1 = .8, 1 = .8, 2 = .5, 3 = .3 1 = 2, 2 = 4, 2 = 7, 2 = 3, 4 = 4. We have deliberately mixed distributions and chosen parameter values such that the transition 0 → 1 → 3 → 5 has much longer waiting time (heavier tailed overall distribution) than the competing path 0 → 2 → 4 → 5.
186
A.V. Huzurbazar / Journal of Statistical Planning and Inference 129 (2005) 181 – 193
The MGF for the overall waiting time using these distributions in (2.2) is 1 1 1 2 − 21 (22 −2s) M(s) = 0.42e −s 1 2 2 1 2 − 21 (22 −2.0s) (2 − s)−1 × 2 1 − 0.5e 2 − s 1 1 1 2 − 21 (22 −2s) −1 + 0.203 e 4 1 − 0.3 1 (1 − s) 1 − s × (3 − s)−1 (4 − s)−1
2 1 2 1 1 2 − 21 (22 −2s) −1 − log 1− 0.3 1(1 − s) − 0.5e 2 − s 1 − s 1 2 2 2 1 2 + 0.15 1 e1 2 − 1 (2 −2s) (1 − s)−1 . 1 − s 2 − s One method for converting MGFs into a densities is to use a saddlepoint approximation. A saddlepoint approximation for the density of one random variable T, is . (2.3) f (t) = {2K
(ˆs )}−1/2 exp{K(ˆs ) − sˆ t}, where K(·) = log M(·), cumulant generating function (CGF) for T , dK , K (·) = ds d 2K K
(·) = and ds 2
K (ˆs ) = t. The last equation, K (ˆs ) = t, is called the saddlepoint equation and sˆ is called the saddlepoint. Using this approximation requires only that the MGF of T exist in an open neighborhood of zero. The literature on saddlepoint approximations in statistics is extensive beginning with Daniels (1954). For an introduction to practical saddlepoint approximations for use in statistics see Huzurbazar, S. (1999). Note that our interest with flowgraph models is to convert the resulting MGF into a density, CDF, etc., and the saddlepoint approximation is one way to do this. There are numerical inversion techniques for Laplace transforms that could also be used. For example, see Ramaswami (2000). We use the saddlepoint with flowgraphs because it is an analytical approximation that is fast, accurate, and well-suited in this case. While expressions for flowgraph MGFs may look complicated, flowgraph algebra lends itself to symbolic computation. The flowgraph structure can be computed using a symbolic algebra package such as MAPLE, which in turn can be used to generate C++ or FORTRAN code for the necessary functions, K(·), K (·), and K
(·), that are inputs to the saddlepoint approximation. The resulting density and hazard functions are given in Fig. 2. The bump in the density is reflective of the manner in which our finite mixture was chosen. This type of modeling allows us to realistically model the shape of the hazard to include portions that are increasing,
A.V. Huzurbazar / Journal of Statistical Planning and Inference 129 (2005) 181 – 193
187
0.08
density
0.06 0.04 0.02 0.00 0
10
20 30 time (in days)
40
50
0
10
20 30 time (in days)
40
50
hazard
0.15 0.10 0.05 0.0
Fig. 2. Density and hazard functions for manufacturing system.
decreasing, or flat (cf. Aalen and Gjessing, 2001). For a given system, there is a wide scope for fine tuning the underlying distributions to reflect the actual system behavior.
3. Construction engineering application A construction engineering project examined the impact of communication and planning on performance (cf. Garcia, 1997). Data on every aspect of operations was collected for twenty firms in New Mexico, USA. Monthly performance is one of the main quantities of interest. Monthly performance is loosely defined as whether the project completed the month on schedule and whether it will complete the next month on schedule. To evaluate completion times, current methods used (reluctantly) by these firms are Critical Path Method (CPM) and Project Evaluation and Review Technique (PERT). Both of these methods are commonly used in management and engineering to identify potential problems in a project plan but it is accepted that they have many shortcomings. See Pritsker et al. (1989) and AbouRizk and Wales (1997) for details. They are used at the beginning of a project and are not very adaptable to the project once the work has begun. They do not provide much assistance to the project manager in terms of timely updates to the project schedule. See Leung and Tam (1999) for details. CPM is also limited to analyzing deterministic systems with no feedback or failures in the project. Projects are always assumed to complete successfully. PERT is limited to using beta and normal distributions. Some current trends include the
188
A.V. Huzurbazar / Journal of Statistical Planning and Inference 129 (2005) 181 – 193
Fig. 3. Flowgraph model for construction engineering project.
use of genetic algorithms (cf. Li et al., 1999). Genetic algorithms and other evolutionary programming methods require large amounts of data for the algorithm’s learning process and a great deal of computation to implement. We use these data to illustrate a Bayesian predictive analysis from the flowgraph model. Fig. 3 shows a flowgraph model for the construction engineering project. State 0 is the planning state and state 1 is the design phase. From state 1, transition is possible to state 6, modifications, or into the monthly planning phase, state 2. The loop from 1 → 6 → 1 represents the process of design with modifications. After the initial planning, the monthly phase begins, represented by the boxed subsystem in Fig. 3. The time frame for the subsystem is one month (=23 working days). State 2 is the monthly planning state, state 3 is the start of the monthly work, state 4 represents delays in the monthly work which can result in modifications to the main project, indicated by a line from 4 → 6, or eventually lead to completion, state 5. Analysis from a flowgraph model gives an entire waiting time distribution as well as the system reliability and hazard functions for any total or partial waiting time. For decision making, this immediately provides the mean, median, or any percentile of the waiting time distribution for a variety of input distributions. The analysis from the current month along with prior information about potential delays gives a predictor distribution for the next month. This gives the project team a figure for how many days behind schedule they are likely to be for the next month allowing them to make necessary modifications such as reorganizing workload or hiring additional workers. Other methods generally only provide the mean and median for only a few distributions. Our data for the twenty firms is restricted to states 2–5 and consist of observations only on the monthly phase so our analysis assumes that from state 4, only transitions to state 5 are possible (p45 = 1). With appropriate data, it would be possible to analyze the entire overall schedule with the flowgraph model given in Fig. 3. Six of the firms experienced no delays. For these firms, we observe transitions from MONTHLY PLAN, to START, to COMPLETION 2 → 3 → 5. Of the 14 firms that experienced delays, 10 were able to recover from these and complete the month on schedule, 4 → 5. The remaining 4 firms were irrecoverably delayed,
A.V. Huzurbazar / Journal of Statistical Planning and Inference 129 (2005) 181 – 193
189
Histogram for 3 > 4
0.12
0.10
0.08
0.06
0.04 Gamma 0.02
0.00
2
6
10
14
18
22
26
30
34
38
Waiting Time Fig. 4. Histogram for data on the 3 → 4.
indicated by the transition 4 → 6. We have no subsequent data for them on the current month and therefore, in the present analysis, we will not analyze this transition further. We use histograms of the data to suggest parametric models to be used with flowgraphs. For example, Fig. 4 presents a histogram of data on the 3 → 4 transition. For illustration, an estimated gamma density is overlaid in the figure. After examining histograms for each state-to-state transition of the data, we make the following parametric assumptions about the branch transmittances: • • • •
Branch 2 → 3, constant 4 days, MGF M23 (s) = exp(4s); Branch 3 → 4, gamma(, ) with mean /, MGF M34 (s) = [/( − s)] , for s < ; Branch 4 → 5, constant 18 days, MGF M45 (s) = exp(18s); Branch 3 → 5, inverse Gaussian (1 , 2 ) with density f (t | 1 , 2 )= √
1 exp{−(21 /t+22 t)/2+1 2 } for t > 0, 1 > 0, 2 0 2t 3/2
and MGF M35 (s) = exp{1 2 −
21 (22 − 2s)},
s<
22 . 2
190
A.V. Huzurbazar / Journal of Statistical Planning and Inference 129 (2005) 181 – 193
Data on 2 → 3 are all within a small fraction of 4 days and data on the 4 → 5 transition are clustered around 18 days. We treat these as roughly constant times to demonstrate how flowgraphs handle constant variables. The 3 → 5 transition has only 6 observations. For illustration, we use an inverse Gaussian distribution to allow for a heavy tail. This is motivated largely by the problem at hand. The 3 → 5 transition means that work is completed without delays. In the construction industry, this is rarely the case. The heavy tail allows for longer waiting times when delays exist but are not recorded. Using (2.1) the MGF of the overall waiting time distribution from 2 → 5. There are two paths through the system, P1 (s) = p34 M23 (s)M34 (s)M45 (s) P2 (s) = p35 M23 (s)M35 (s) where p34 + p35 = 1, and we take p12 = 1 and p45 = 1 for the monthly phase analysis. The MGF of the total waiting time distribution for monthly time to completion is M25 (s) = e p34 4s
−s
e
18s
2 2 + e (1 − p34 ) exp 1 2 − 1 (2 − 2s) (3.1) 4s
for s < min(, 22 /2), 1 > 0, 2 > 0, > 0, > 0, 0 p34 leq1. The model can accommodate any transition time distribution, discrete or continuous. In particular, if 4 → 5 took on the values 18 and 19 with probabilities 23 and 13 , the term e18s in (3.1) would be replaced by (2/3)e18s + (1/3)e19s . Our goal is to compute the Bayes predictive density of the overall waiting time from planning to completion, 2 → 5, for a given month. We let T represent this random time. The Bayes predictive density of a future observable Z given data D about a system with parameter vector is f (z | D) =
f (z | )L( | D)() d ≡ E( | D) {f (z | )}, L( | D)() d
(3.2)
where z has density f (z | ), L( | D) is the likelihood function, and () is the prior. Our parameter vector is = (p34 , , , 1 , 2 ). To compute the Bayes Predictive density of (3.6), we assign a U nif (0, 1) prior for p34 and we assume diffuse U nif (0, 20) priors on all of the other parameters. The likelihood is separable in p34 , (, ), and (1 , 2 ). This results in a beta(15, 7) posterior distribution on p34 and the posterior on the remaining parameters is given by
(, , 1 , 2 | D) ∝
()
14 14
ti−1
e−
14
i=1 ti
i=1 20 20 2 2 1 × 61 exp − 1 − 2 t j + 6 1 2 2 tj 2 j =15
j =15
(3.3)
A.V. Huzurbazar / Journal of Statistical Planning and Inference 129 (2005) 181 – 193
191
0.04
density
0.03
0.02
0.01
0.00 0
10
20
30
40
50
t
Fig. 5. Bayes predictive density function for the completion time of the construction engineering project.
for 0 < , , 1 , 2 < 20. where the ti , . . . , t14 are the observed waiting times for the transition 3 → 4 and t15 , . . . , t20 are the observed waiting times for the transition 3 → 5. Note that the waiting times for transitions 2 → 3 and 4 → 5 are constant. Monte Carlo methods are used to perform the integration in (3.2) using a random sample of values from the posterior in (3.3) and from the beta(15, 7) posterior. Posterior sampling is done via slice sampling (cf. Neal, 2003), a Markov chain Monte Carlo method. Computationally, this involves using slice sampling to get a random sample 1 , . . . , m from the posterior. For each k value, we evaluate fˆ(z | k ) and L(k | D). fˆ(z | k ) is computed using the saddlepoint approximation (2.3) to convert the underlying flowgraph model MGF (3.1) into a density value. Given data D, the likelihood function for each k is evaluated based on the parametric assumptions on the branch waiting times. Figs. 5 and 6 show the Bayes predictive density and hazard function for this project. The predictive mean is 32.29 days and the standard deviation is 141.91 days so we do not expect to remain on schedule. The large standard deviation is reflective of the very heavy tailed distribution for the 3 → 5 transition and the fact that the once the project falls behind schedule it takes much longer to complete due to compounding problems. The hazard increases steadily until about 30 days and then drops off. This is reflective of the fact that project delays make it hard to complete in the 20 additional days shown in the figure. However, once the problems are resolved, the the hazard for the completion time will start to increase again in the furture. The 25th, 50th, 75th, 90th, 95th, and 99th percentiles are: 20.99, 27.28, 34.25, 44.57, 55.33 and 103.76 days. Any other desired quantity is also easily obtained from the predictive density. Alternatively, we could have used maximum likelihood to estimate the density of Z as f (z | ˆ ). The flowgraph provides the MGF which can be evaluated at ˆ , the MLE, and inverted to give an approximate density, fˆ(z | ˆ ), using the saddlepoint approximation. The
192
A.V. Huzurbazar / Journal of Statistical Planning and Inference 129 (2005) 181 – 193
0.10
hazard
0.08
0.06
0.04
0.02
0.00 0
10
20
30
40
50
t
Fig. 6. Bayes predictive hazard function for the completion time of the construction engineering project.
shape of the hazard is increasing and decreasing at varying rates. This is typical of hazards from flowgraph models since we are modeling the actual waiting times and not making assumptions about the hazards.
4. Conclusion Flowgraphs provide a data-analytic approach for modeling time-to-event data. Computational aspects begin with determining a flowgraph model for a given system. Estimation based on the flowgraph model can be performed in the Bayesian framework or via maximum likelihood. MGFs of total waiting times of interest are then computed from the flowgraph model. Inversion methods such as saddlepoint approximations can be used to convert this into a density, CDF, reliability, and hazard functions. Flowgraph models allow data analysis for semi-Markov processes using non-exponential waiting times. They also allow analysis of systems with feedforward and feedback loops in the presence of general waiting time distributions. Flowgraphs model semi-Markov processes that are applicable to a wide variety of situations. In the engineering setting, they can be used for systems analysis and design of complex systems.
Acknowledgements The author thanks Sandra Garcia, Federal Highway Administration, Department of Transportation, Boise, ID, for providing the construction engineering data and two anonymous referees for helpful comments.
A.V. Huzurbazar / Journal of Statistical Planning and Inference 129 (2005) 181 – 193
193
References Aalen, O., Gjessing, H., 2001. Understanding the shape of the hazard rate: a process point of view. Statist. Sci. 16, 1–13. AbouRizk, S., Wales, R., 1997. Combined discrete-event/continuous simulation for project planning. J. Construct. Eng. Manage. 123, 11–20. D’Azzo, J., Houpis, C., 1981. Linear Control System Analysis and Design: Conventional and Modern, McGrawHill, New York. Butler, R., Huzurbazar, A., 1997. Stochastic network models for survival analysis. J. Amer. Statist. Assoc. 92, 246–257. Daniels, H., 1954. Saddlepoint approximations in statistics. Ann. Math. Statist. 25, 631–650. Dorf, R., Bishop, R., 1995. Modern Control Systems, Addison-Wesley, Reading, MA. Gajic, Z., Lelic, M., 1996. Modern Control Systems Engineering, Prentice-Hall, New York. Garcia, S., 1997. Planning and Communication in Construction: Impacts on Performance. M.S. Thesis, University of New Mexico. Huzurbazar, A., 1999. Flowgraph models for generalized phase type distributions with non-exponential waiting times. Scand. J. Statist. 26, 145–157. Huzurbazar, A., 2000. Modeling and analysis of engineering systems data using flowgraph models. Technometrics 42, 300–306. Huzurbazar, S., 1999. Practical saddlepoint approximations. Amer. Statist. 53, 225–232. Leung, A., Tam, C., 1999. Models for assessing hoisting times of tower cranes. J. Construct. Eng. Manage. 125, 385–391. Li, H., Cao, J., Love, P., 1999. Using machine learning and GA to solve time-cost trade-off problems. J. Construct. Eng. Manage. 125, 347–353. Mason, S.J., 1953. Feedback theory—some properties of signal flow graphs. Proc. Inst. Radio Eng. 41, 1144–1156. Neal, R., 2003. Slice Sampling. Ann. Statist. 31, 705–767. Ouhbi, B., Limnios, N., 1999. Nonparametric estimation for semi-Markov processes based on its hazard rate functions. Statist. Infer. Stoch. Process. 2, 151–173. Ouhbi, B., Limnios, N., 2003. Nonparametric reliability estimation of semi-Markov processes. J. Statist. Plann. Infer. 109, 155–165. Pritsker, A., Sigal, C., Hammesfahr, R., 1989. SLAM II Network Models for Decision Support, Prentice-Hall, New Jersey. Ramaswami, V., 2000. Algorithmic analysis of stochastic models: the changing face of mathematics. Ramanujan Endowment Lecture at Anna University, Chennai, India. Whitehouse, G., 1973. Systems Analysis and Design Using Network Techniques. Prentice Hall, NJ.