(navigation image)
Home American Libraries | Canadian Libraries | Universal Library | Community Texts | Project Gutenberg | Children's Library | Biodiversity Heritage Library | Additional Collections
Search: Advanced Search
Anonymous User (login or join us)
Upload
See other formats

Full text of "Stochastic modeling : ideas and techniques"

LIBKARY 

RESEARCH REPORTS DIV 
NAVAL POSTGRADUATE 
MONTEREY, CALIFORNIA 



NPS55-82-005 






NAVAL POSTGRADUATE SCHOOL 

Monterey, California 




STOCHASTIC MODELING: IDEAS AND 
TECHNIQUES 

by 
Donald P. Gaver 



January 1982 



Approved for public release; distribution unlimited 

Prepared for: 

PK ^°f of Naval Research 

ngton, Virginia 22217 



FEDDOCS 

D 208.14/2:NPS-55-82-005 



NAVAL POSTGRADUATE SCHOOL 
MONTEREY, CALIFORNIA 

Rear Admiral J. J. Ekelund David A. Schrady 

Superintendent Acting Provost 

Work on this report was sponsored in part by the Office of Naval Research, 
Arlington, VA. 

Reproduction of all or part of this report is authorized. 



Unclassified 



SECURITY CLASSIFICATION OF THIS PAGE (When Data Entered) 



REPORT DOCUMENTATION PAGE 



READ INSTRUCTIONS 
BEFORE COMPLETING FORM 



I. REPORT NUMBER 

NPS55-82-005 



2. GOVT ACCESSION NO 



3. .RECIPIENT'S CATALOG NUMBER 



4. TITLE (and Subtitle) 

STOCHASTIC MODELING: IDEAS AND 
TECHNIQUES 



5. TYPE OF REPORT 6 PERIOD COVERED 

Technical 



6. PERFORMING ORG. REPORT NUMBER 



7. AUTHORf*; 

Donald P. Gaver 



8. CONTRACT OR GRANT NUMBERf*,! 



9. PERFORMING ORGANIZATION NAME AND ADDRESS 

Naval Postgraduate School 
Monterey, CA 93940 
University Libre de Bruxelles 
Brussels, Belgium 



10. PROGRAM ELEMENT. PROJECT. TASK 
AREA 4 WORK UNIT NUMBERS 

61153N; RR14-05-OE 
N000 1482WR20017 



II. CONTROLLING OFFICE NAME AND ADDRESS 

Naval Postgraduate School 
Monterey, CA 93940 



1? <EPORT DATE 

January 1982 



13. NUMBER OF PAGES 
82 



U. MONITORING AGENCY NAME a ADORESSf// dl Iterant Irom Controlling Oltice) 

Chief of Naval Research 
Arlington, Virginia 22217 



15. SECURITY CLASS, (ol thle report) 

Unclassified 



15a. DECLASSIFICATION/ DOWN GRADING 
SCHEDULE 



16. DISTRIBUTION STATEMENT (ol thia Report) 

Approved for public release; distribution unlimited 



17. DISTRIBUTION STATEMENT (of the abstract entered In Block 20, It dltlerent Irom Report) 



18. SUPPLEMENTARY NOTES 



19. KEY WORDS (Continue on reverse aide It neceaaary and Identity by block number) 



applied probability 
probability modeling 
reliability 
availability 
waiting lines 



20. ABSTRACT (Continue on reverae aide It neceaaary and Identity by block number) 



This report summarizes the contents of lectures given on prob- 
ability modeling and reports some new results on the availability 
of inspected systems of redundant systems in random environments, 
and on "sculptured distributions". 



dd ,; 



FORM 
AN 73 



1473 EDITION OF 1 NOV 65 IS OBSOLETE 

S/N 0102-014' 6601 | 



Unclassified 



SECURITY CLASSIFICATION OF THIS PAGE (When Data Bntered) 



Stochastic Modeling: Ideas and Techniques 
Donald P. Gaver 

1 . Introduction 

The primary purpose of this chapter is to summarize the 
contents of lectures on stochastic modeling presented at the 
Universite Libre de Bruxelles (ULB) in the period March-May, 
1981. Much of the material selected for presentation was from 
the standard menu of probabilistic topics typical of a second 
course as given to engineers, operations researchers, 
statisticians, or computer scientists. An attempt was made to 
emphasize a modeling attitude rather than details of mathematical 
rigor, illustrating with problems and techniques that are not 
often prominent in such courses. For example, attention was 
given to problems of, and models for, redundant system reliability 
and availability, queueing with priorities, first-passage times 
and areas under path functions of stochastic processes, (total 
waiting times), and various other topics. Also included was a 
brief account of aspects of modern data analysis, with the 

implication that its usefulness is significant at the pre-modeling 
and model-assessment stage of an investigation. 

A secondary, but gratifying, purpose is to briefly report 
on cooperative work initiated with faculty and students at ULB. 
I wish to mention the enjoyable collaboration with Dr. Guy Latouche 
on development of efficient computational methods for repairman- 
like Markov models in random environments, and with Ph. Collard 
on the application of sculptured distributions in the simulation 



evaluation of certain scheduling algorithms. The interest and 
warm hospitality of Prof. Guy Louchard, head of the Dept.Of 
Computer Science at ULB, was also much appreciated. 



2 . The Total Modeling Process: Brief Overview 

It is coming to be recognized that the topic of mathematical 
modeling (including stochastic modeling) exists in its own right 
as a subject suitable for a formal university course; see 
Bender [1978], The modeling step is part of a process of several 
stages or steps; these may be expressed as follows (Gaver and 
Thompson [1973] ) : 

(a) Identify the general problem area or s ituatio n; identify 
specific questions concerning that area. 

(b) Obtain and analyze subject-matter information and data 
relating to the problem area. Often an examination 

of such information and data will suggest suitably 
formulated questions, as in (a), 
(c). Construct a preliminary model, or models, representing 
the important features of the situation. Deduce some 
model implications. 

(d) Refer the result of (c) to subject-matter specialists 
and decision-makers for qualitative critique; revise 
the model accordingly. This likely means re-doing 

(a) - (d) . 

(e) Assess the empirical validity of the model to the degree 
possible. Check the sensitivity of model conclusions 

to changes in model assumptions (sub-model inputs) , 
and to data variations. Submit to judgement by 
subject-matter experts -- but anticipate differences 
of opinion! The modeling, and re-modeling, process 
may help to reconcile such differences. 



(f) Compute required answers to interesting questions. 
Assess the degree of uncertainty in these answers 
possibly resulting from model mis-specification, data 
bias or other deficiency, computational error, and 
sampling error in estimates of basic parameters or 

in simulation results used to supply model implications, 

(g) Communicate, and aid in implementing, the results of 
the model. 

(h) Monitor the situation for possible changes in the 
environment, and hence for the necessity to change 
the model. 
Of course the emphasis in these notes (and in the lectures, 
was) upon the actual modeling step, (c) . However, some attention 
was given to the display of data for pre-model examination 
(Tukey's exploratory data analysis), and to model parameter 
estimation techniques, particularly those robust methods that 
attempt to deal with questions of data deficiencies. 



3 . Topics in Outline 

In this section we out-line the basic contents of the 
lectures. These were in general arranged so as to first present 
mathematical definitions and properties, and then to illustrate 
in terms of sample models for various situations. 
( 1 ) Review of Probabilistic Concepts, Particularly Conditioning . 

In this lecture the following basic notions of probability 
were defined or reviewed: random experiment or trial, sample 
or event space, events and combinations of events, probability as 
a function with rules for combination, conditional probability, 
independence, and Bayes 1 Theorem, random variables and their 
moments or expectations, transforms (characteristic function, 
Laplace transform, and generating function) and their moment- 
generating, convolution, unicity, and continuity properties, plus 
properties of conditional expectations. In addition, certain 
classical univariate distributions were reviewed (Normal/Gaussian, 
log-Normal, exponential and gamma, etc.) 

By way of illustration, a simple problem of equipment (or 
possibly software) unrealiability was considered. 

Situation: Suppose a system is made up of components that 
individually fail after a time because of the action of faults ; 
the latter may be the result of component mis-design, or attributable 
to bad installation or adjustment ("human error"), or to a mistake 
in computer program coding. We wish to relate system failure rate 
to initial fault content. 

Model: N is a random variable (rv) representing the 
number of faults initially installed unwittingly in the system. 
Let {T., i = 1,2, ..., N} be the sequence of rv describing 



the failure time of each fault, measured from time at which the 
system begins use; T. may actually be the time at which the 
service of the particular component is first requested. Suppose 
the system fails at 

X N = min {T-, T 2 , ..., T„}. (1.1) 

Under very simple conditions, namely that all components have 
the same distribution, F(t), of failure time, and all failure 
times are independent, simple conditional probability arguments 
yield 

P(X N > t} = E N [(l-F(t)) N ] (1.2) 

where g is the generating function of the number of faults 
originally sown in the system, and F (t) = 1-F (t) is the 
survival time distribution, per fault. It is easy to see that 
P{X = 00} = p{n = 0}, possibly > , so the derived distribution 
of X is quite possibly dishonest. Note that while in general 
explicit expressions for expectations cannot be obtained (may 
not even exist), such summaries as the median, 90% point, etc., 
may if g (F(t)) can be explicitly inverted, e.g. for 
N ~ Poisson and F ~ Exponential. The simplistic assumption 
of the model may be relaxed, allowing for different T. distri- 
butions, dependence, and so on, and an additional random death 
time, D, applying to the total system can be introduced to in- 
duce eventual failure (of physical equipment), or biological 
death in finite time. There will be less analytical tractability, 



but simulation may be used to assess system behavior. Statistical 
estimation problems may be addressed as well; a suitable version 
of (1.2) will provide a likelihood function. 

Another example of the applicability of a simple conditioning 
argument is the following. 

Situation ; When an individual speaks on a telephone or 
telecommunication channel the conversation is an alternating 
sequence of talk-spurts and pauses. Similarly, a job being 
processed on a computer goes through an alternating sequence of 
CPU (compute) times and 10 (input-output) times. Model the 
total time of the conversation or job processing time, and 
particularly the joint distribution of busy and idle segments. 
Model 1 : Let {X., i = 1,2, ..., K} and {Y. , i = 1,2, ..., K) 
be the durations of talk spurts and pauses, respectively, and 
let K be the number of each. The simplest model assumes 
{X.} and {Y.} to be independently and identically distributed 
(IID) sequences of rv , and themselves to be conditionally in- 
dependent, given K, also a rv. The joint distribution of total 
talking (or processing) time, X, and total pause time (10 time), 
Y, is thus, by simple conditioning, 

? k* k* 

P{X < x, Y < y} = I F v (x) • F v (y) K P{K=k} (1.3) 

k=l X Y 

The joint Laplace transform is (s. , s_, > 0) 



— s X — s Y 

e 2 j I [F X ( S;L )F (s 2 )T P(K=k) 



1" "2" r r ~ , , ~ , ,ik 



k=l 



(1.4) 



= %^X (S 1^Y (S 2 )] ' 



where g v is the generating function of K. Put s n = s„ = s 
is. 12 

to recover the transform of X+Y = L, the total conversation 
length. In case 

X. - Expon(A) and Y. Expon(y) and K Geom(a), independent 



and 



" S 1 X " S 2 Y 
e e 



- I 
k=l 



A+s 



V 



lP+s J 



(1-a) a 



k-1 



Ay (1-a) 



(A+s ) (y+s ) - Aya 



(1.5) 



E[X] = [A (1-a)] 1 , E[Y] = [y(l-a)] 1 



(1.6) 



Furthermore (a = 1-a) .' 



E[L] = E[X] + e[y] = 



A+y 
Aya 



(1.7) 



and 



Var[L] = (E[Lj) 



Aya 



(1.8) 



Notice also that the mechanism of randomization of a sum , or 
mixing (see Feller [1966])which has given (1.3) may be used to 
generate families of bivariate (multivariate) exponential 
distributions for other modeling purposes. 

Model 2 : A plausible alternative to the above model assumes 
X. and Y. are not independent, being possibly positively 
correlated -- a long talkspurt tending to result in a long 
pause (response by conversationalist). Most simply, 



Y. = $X . , 3 > 0. Then again transform in the X - Expon(A) case 
to get 

_<= Y -c v" 

A (1-a) 



- S 1 X " S 2 Y 
e e 



s x + 6s 2 + A (1-a) 



or if 3 = A/p which preserves the marginals of Model 1, 



A (1-a) 



s, + — s„ + A (1-a) 
1 y 2 



Aya 



s,y + s A + Aya 



(1.9) 



It is now immediate that the marginal distribution of 

X ~ Expon(Aa) , Y - Expon(ya) and that now the df of the 

Aya 



total time, L, is simply Expon 



[A+y 



-- a much simpler form 



than that occurring in Model 1 above, which involves a Bessel 
function. The variance of L in Model 2, being (e[l]) , is 
also larger than that for Model 1, see (1.8), suggesting that 
the former model has a longer tail, hence predicting a greater 
proportion of extremely long conversations. 

The above illustrates that the same situation can easily 
give rise to two — or more — different models, depending 
upon the manner in which stochastic assumptions are introduced. 
At best, the introduction should be guided by observed data; 
at least, sensitivity analyses using different assumptions can 
outline the range of specification uncertainty. 



(2) Models Involving Repeated Trials 

A great many situations may be initially modeled in terms of 
repeated independent trials, where this means that on each of a 
possibly countably infinite number of occasions a trial (or 
experiment, or observation) is performed, with outcome X. 

(possibly a vector random variable (rv) ) on the ith trial; 
{X., i = 1, 2, ...} are IID rv. Bernoulli trials are a prime 
example: flip a biased coin indefinitely; let I. be one if a 
Head (Success) results, and zero if a Tail (Failure) otherwise, 
and assume that probability of success on any trial is independent 
of all previous outcomes. Equally, X. may be the winnings on a 
bet at occasion i, with X. in dollars and either positive or 
negative. Or X. may even represent the increase or decrease 
in a common stock price on the New York Stock exchange, according 
to some observers. 

It is convenient, but somewhat more questionable, to adopt 
the repeated trial model for modeling real operational and 
physical phenomena, yet it is often done uncritically. 

For instance the lifetimes (times between failure) of 
computing equipment are often modeled by IID rv, repair times 
likewise, queueing system inter-arrivals and service times as 
well, inventory demand sizes, sizes of deposits of resources 

(petroleum) as well, ... the list is very long, and observational 
support for these assumptions is usually conspicuously lacking. The 
Attraction of the repeated trials model is mainly its mathematical 
tractability , which leads to elegant and appealing results. 



10 



Often a brief data analysis in terms of marginal distributions 
of observed X. seems to provide justification. The simple 
repeated trials model cannot well represent, say, systematic 
daily changes in job inter-arrival times, or numbers of jobs 
per hour, at a computer center, seasonal effects on such computer 
system demand, or the influence of other variables such as the 
introduction of a new class of computer users upon a measure of 
computer loading. Also, the model does not well describe the 
sequence of daily rainfalls in a region, nor many other environ- 
mental variables. Some examples follow in which repeated trial 
models seem initially plausible, but which doubtless can stand 
improvement. 

Situation ; A structure is to be designed to withstand (wind- 
generated, or seismic or other) shocks. Question: how long will 
it survive the environment stresses, given its initial strength Z? 
Model: Let the structure have strength Z (suitable units). If 
the ith year's maximum shock is X., assume {X.} to be IID with 

df F (x) . Then the time T until structure failure exceeds t 
x 

(t = 1, 2, ...) iff VX ± < Z, i = 1,2, ..., t, so 

P{T > t|z} = [F x (Z)] t , (2.1) 

the geometric distribution, 

while if Z itself is regarded as random 



P{T > t} = E z P{T > t|z) = E z ([F x (Z)] t } 



(2.2) 



^ (E z {[F x (Z)]}) t ; 



11 



the unconditional distribution is a convex (probability) combin- 
ation of exponentials. It will resemble an exponential, but often 

has an extended right tail. It is obviously important that the 
condition on equipment stress, Z, be removed at the appropriate 
stage -- at the end. Removal of the condition before each 
"strength test" would be appropriate only if the structure were 
completely repaired or replaced with another having the unchanged 
distribution of Z before each trial. 

Note that assuming yearly maximum environmental events to 
be IID is intuitively plausible, but some statistical evidence 
exists for truly long-range correlation in weather data that 
may call even this assumption into question. 

The Bernoulli counting process is an important special case 
of repeated trials, on each of which either Success (probability 
p) or Failure (probability q = 1-p) occurs; {N(t), t = 1,2, ...} 
is the number of Successes in d,t]). Times (number of trials) 

t between successive Successes are IID and geometrically distrib- 

k-1 
uted (P{t = k} = q p) . The number N(t) of Successes in 

fixed time is Binomial (p , t) , and is in turn approxi- 
mately Normal (tp, tpq) as t •+ °°. Of course Bernoulli trials 
describe the outcomes of many other repeated trial situations: 
for instance, the number of jobs submitted to a processing 
facility requiring more than x time units of processing may 
be modeled as a Bernoulli counting process with p = 1 - F (x) , 
F being the distribution of job processing time. 

Generalizations of Bernoulli trial situations may be (a) 

to variations of success probability with trial number ("time", 

t 
or "space"): P{X. = 1} = p.. Here E[N(t)] = I P- and 

11 i=l 1 

t _ ! t 

Var[N(t)] = I p.q. £ tp(l-p), where P = 7" I Pi' so variability 
i=l i=l 

is u nder-represented if trial-to-trial probabilities change, but 

12 



a Normal approximation to N(t) may still hold. A second qeneralization 
is to (b) independently randomize p in the Binomial distribution, 
say according to a Beta distribution, creating the Beta-Binomial dis- 
tribution. This has found useful in reliability modeling and 
in Bayesian inference. A third and important elaboration, (c) , is to 
develop a statistical regression model for success probability p 
based conveniently on the logist ic function: 

a+Bu. 

P- = — T5 — ; (2.3) 

*i a+3u. 

1+e 1 

here u. (possibly a vector) represents the influence of other 

factors upon success probability. Given observations of the 

form (I., u.), where I. = 1 indicates success on trial i, one 
11 l 

can estimate a and 3 (vector) by maximum likelihood; see 
Cox [1969] • Generalizations to multiple-category situations are 
possible, and computational methods for parameter estimation 
and model assessment have been devised; see Pregibon [1901J. 

The distribution of the number of counts N(t) in t trials 
of a Bernoulli trials process can be computed by making use of a 
forward equation . Let 

P. (t) = P{N(.t) = j |N(0) = 0}. 

Then 

P.(t) = P. (t-1) • (1-p) + P._ 1 (t-1) -p;(0 < j It); (2.4) 

on the basis of conditioning on events that have happened up to 
t-1. One can generalize to non-stationary success probabilities 
easily: 



13 



P. (t) = P. (t-1)- (1-p ) + P. (t-l)p. . (2.5) 

J J T. J X XL 

Initial conditions may be 

P Q (0) = 1, P. (0) =0, j = 1,2, ... 

A further generalization allows success probability to depend 
upon the number of previous successes; then 

P j( t) = P j( t-l)U-p. <t ) ♦ P^U-Up^^ , (2.6) 

and the distribution P. (t) can easily be computed recursively 
given the success probabilities 

p = P{N(t) = j+l|N(t) = j} . (2.7) 

This is a preview of ideas of Markov chains, to be treated later. 
These expressions are introduced to suggest early that the answers 
to interesting and comparatively complex problems can be directly 
computed numerically (in this case iteratively, starting from 
the initial conditions) . Closed- form expressions such as the 
Binomial distribution are handy, and Normal approximations are 
even handier, but one need not modify the facts merely for the sake 
of convenience. 



14 



(3) Sums of Repeated Trial (IIP) Random Variables; "Large 
Deviations " 

Models for total demand for physical inventory or for 
facility (computer) time often naturally involve sums of varying 
components, modelled as rv. ; thus total demand from n sources, 
or over n time periods, is 



S = X, + X n + . . . + X . 

n 1 2 n 



(3.1) 



If X^ is the (dollar) profit in the i th year for some enterprise, 
then a financial measure of success is 



n 

S (r) = y X.r 1 
n ■ -, i 
i=l 



(3.2) 



where r is a discount rate (0 < r < 1) 



Situation . A computer center experiences varying monthly demands, 
X. for the ith month. Here are answers to several simple 
questions involving sums of X.s. 

The expected yearly (n-period) demand is 

n 



E[S ] = I E[X. ], 
n- 1 , L - l 

i = l 



(3.3) 



the sum of the expected monthly demands, and also 

n 

Var[S ] = I Var[X.] 
n . , i 
i=l 



(3.4) 



provided the X.s are uncorrelated . Importantly, as n 



F ,(x) 

n 



P< 



S - E[S ] 
n L n 

/VarLS J 
n 



< x 



> % / 



1 2 
"2 Z 



= $ (x) , 



'2n 



(3.5) 



15 



i.e. S becomes approximately Normally distributed no matter 
what the distributions of X . , by the Central Limit Theorem, 
provided the X. components are all of about the same size 
(certainly if they all come independently from the same parent 
distribution with finite mean and variance). For smallish n or 
distinctly non-Normal components the approximation is improved 
by an Edaeworth expansion (Feller [1966]), wherein for the eaual 
component example 



F gI (x) = $(x) + 
n 



ry 3 



3 



(x 2 -l) |* + R (3.6) 

dx n 



where R = 
n 



1 



/n~ 



densities. Here y = E 



6/n 

, and the components are assumed to have 

3 



(X - E[X]) 



y 3 
and the term — ^ = y 



G 

is the conventional dimensionless skewness measure for a dis- 
tribution, being zero for symmetric distributions (Normal) , 
and being +2 for the Exponential. Additional terms involving 
kurtosis (4th moments) improve the approximation, but it is 
possible that Edgeworth numerical values can be "infeasible" : 
the approximation can actually decrease with x in certain 
ranges. Nevertheless the Edgeworth series has been usefully applied, 
even to unequal component situations, for estimating the loss 

of capacity of an electric utility; see Levy & Kahn [1981]. 

A useful alternative is the method of large deviations , 
Feller [1966], Daniels [1954], and others. The ingenious 
idea is to tilt (or sculpture) the df. components so as to make 
a Normal approximation more effective at predicting the prob- 
ability that S > x for large x. For equal components with 



16 



d.f. F(x) look at the tilted probability measure (assumed to 
exist for s > 0, which sometimes restricts the theoretical ap- 
plicability) : 

v(dx) = eSXF(dx) n 7> 

viaxj , (3.7) 

e 

/v sX 

ip(s) = £n F(s) = in E[e ] being a cumulant generating 

function for F, or X. Manipulations show that 



'{S n > z} = / F n *(dz) = e n ^ (s) J e" SX V n * (dx) , 



(3.8) 



n* 
and the idea is to approximate V by a Normal centered at z, 

a feat that can be accomplished by choice of s. It turns out 

that it is necessary to solve (sometimes numerically) for s(z) 

the equation 

z = nip' (s) (3.9) 

in order that the mean of the approximating Normal be at z ; 
the variance is nij;"(s). Finally, 



P{S n > z} % exp n{i|>[s(z)] - s(z)^'[s(z)] + | s 2 (z) ty" [ s (z ) ] } x 

1 2 

oo -—v 



i- / e 2 dv . (3.10) 



s(z)/mj/'Ls(z) j 

The above technique can also be applied to a compound 
Poisson model (n is replaced conditionally by N(t), the 
counting process of a Poisson process, and the condition then 
removed) . Such models are frequently employed in inventory 
studies; apparently the large deviations approximation has not 
been applied in that area. 

17 



( 4 ) Bernoulli Trials and Poisson Process: Rare Events 

Bernoulli Trials are a special case of the Repeated Trials 
model, with events occurring ("Success") or not permitted to 
occur ("Failure") at specific integer time points, often equally 
spaced. In practice the fixed intervals between trials may 
be largely arbitrary, and it is attractive to think of events 

occurring at any (real-valued) time; from this comes the Poisson 
process. One approach to the P.P. properties is to consider a 
B.T. process to operate over time t with unit time steps, and 
then refine the time steps (e.g. let t = 1 day and starting with 
possible demands at 15-minute intervals, then down to 7.5 minutes, 
then to 3.75 etc. ) to create a sequence of B.T. models. The 
limit of the sequence after ultimate refinement describes the 
P.P. 

Specifically, let T(k) be the generic time between successes 
in the (kth) B.T. model with time steps 1/2 (k = 0,1,2, ...). 
This means that the actual number of steps in time t for B.T. 
model k is 2 t; correspondingly, let the probability of success 
per step be p/2 . By conditioning on the first step's outcome 
this means that 

E[T(k)] = l/2 k + (1 - p/2 k ) • E[T(k)] (4.1) 

so E[T(k)] = — for every model, as should be true. Furthermore, 
as k -> °° so time steps become arbitrarily small, 

P{T(k) > t} = (1 - p/2 k ) t ' 2 ■* e" tp (4-2) 

inter-event times become exponentially distributed in the (P.P.) 



18 



limit. Furthermore the number of P.P. events ("Successes") 
in time t have the Poisson distribution. 

The P.P. is usefully invoked for many modelling purposes. 
Situation . Consider a sequence of days on which demands for 
computer service (time) are made, and focus on the occurrence 
patterns of runs (uninterrupted sequences) of h igh-de mand days. 
Question: what is the distribution of times between successive 
runs, and what is the distribution of the number of such runs 
in a fixed time t? It will turn out that if either the run 
lengths are long, or if the probability of a high-demand day 
is small, that runs tend to occur as a Poisson process if the 
time scale is appropriate. 

Model . Begin by modelling individual high demand day occurrences 
as successes in B.T. Let t, (k) represent the time until the 
1st occurrence of a run of length k, and, measured from the end 
of such a run, let t, (k). , t-UOi ... t^Oc),.. be the time until 
the 2nd, 3rd, ... ith, such run is realized. By the B.T. as- 
sumption {t.. (JO, i = 1,2, ...} is an IID sequence of rv. 
Then we can represent x (JO by conditioning on the events that 
may occur in the first k trials: 

with prob. p 



x(k) 



1 + t' (k) with prob, 



j + t 1 (k) with prob. p- 1 q 



k-1 
k + t' (k) with prob. p q 



(4.3) 



19 



where x 1 (k) is an independent replica of any x (k) : the idea is 
that the process starts over once a failure occurs to spoil a run. 
Alternatively , 



T(k) = \ 



k with prob. p 

R(k) + t' (k) with prob. l-p ] 



where 



j-l 
P{R(k) = j} = 23^— , 
1-p 



3 -Lf '/ •••/ K , 



(4.4) 



(4.5) 



a truncated geometric distribution. From these come the 
generating function of t (k) , and in principle its distribution: 



[z T < k >] . 




k k 
z p 


1 - 


- qz 


r i-(pz) k l 

l-pz 



Differentiation gives the mean 



E[x(k)] = k + ^ 



P 



1_ kp^ 

1-p , k 

1-p _ 



Qj 



P K (1"P) 



(4.6) 



(4.7) 



the approximation holding if either p -> or k ■* °°. In either 
case the run is a rare event. 

While explicit inversion of the expression for E[z ] 
is possible by use of partial fractions, the result is quite 
complicated. On the other hand, look for the distribution of 



T*(k) = t (k)/E[x (k)] 
when E[x(k)] becomes large. The expectation 



(4.8) 



20 



E[e- ST * (k) ] = E 



e -ST (k)/E[T(k) ] 



(4.9) 



can be obtained from the generating function (by putting 
z = exp[-sp q]); next let either p-*0 (rare individual events) 
or k->-°° (long runs) to find that this transform converges to 
(1+s) . Then by the unicity theorem for transforms (Feller [1966]) 
the normalized rv t* (k) is approximately unit exponentially 
distributed, i.e. 

P{x(k) < t E[-r(k)]} % 1 - e _t (4.10) 

and furthermore the distribution of the number of k-runs in time t, 
N, (t) , is approximately Poisson. Deviation from the Poisson 
(indicated by over-variance) may signify that the underlying 
demand generating process is inhomogeneous or cluster-prone in 
time, and that extra facilities are required to reduce backlogs. 
Examination of runs is one way to check the validity of the 
basic modeling assumption of Bernoulli trials. 

Similar limiting arguments simplify other situations in- 
volving rare events that are generated by even more complicated 
processes. See work on first-passage times for combinations of 
random loads by Gaver, Jacobs and Latonche [1981], 



21 



(5) Markov Models: General Comments 

The basic theory of Markov chains and processes, both in 
discrete and continuous time, is well introduced in standard 
texts such as Feller [1966], Chung [1967], Karlin & Taylor [1975] 
Kleinrockf 1976], and needs no systematic coverage, only review 
and illustration. By way of review, recollect the ideas of 
various possible state space definitions: integers, integer 
and real numbers ( "ages" ) , real numbers (e.g. virtual waiting times 
in queues) ; times (index sets) either discrete and 
equally-spaced or imbedded or continuous time; Markov property 
defined by conditional probabilities ("The future is independent 
of the past, given the present"). Carry on to matrix represen- 
tation of the state probabilities after t (0,1,2,...) time 
steps, forward and backward Chapman-Kolmogorov equations, 
generalize to discrete state Markov chain in continuous time 
with exponential sojourns in states, state classification 
emphasizing irreducible chains and transient chains (with at 
least one absorbing barrier) , recurrent events and first-passage 
times and absorption probabilities, generating functions and 
other transforms. 

Simple Markovian assumptions, i.e. that a scalar state 
rv X(t)f where t is time or space, is Markov, introduce de- 
pendence in a plausible and tractable manner. Usually it is 
necessary to assume, for example, that the one-step transition 
probabilities (discrete state, discrete time): 

p ± . = P{X(t) = j |x(t-l) = i}, (5.1) 



22 



are time-homogeneous in order to obtain explicit neat solutions. 
Analogous assumptions must be made about discrete-state Markov 
processes in continuous time, wherein A. is the rate of de- 
parture from state i (exponentially distributed sojourn time 
parameter), and p. . is the corresponding probability of move 
from i to destination j. Of course a known deterministic 
time dependence, involving daily or weekly cycles, and trends 
can be dealt with by numerically multiplying the transition 
probability matrices. 

More irregular changes in process behaviour can be repre- 
sented as the effect of randomly changing external events, or 
rando m environments for short. In such models the actual 

primary process transition parameters (e.g. p. ., or A.) 

-*- J 

change in time under the influence of such environmental 
factors as seismic vibration, temperature and humidity, ocean 
sea state, wind speed or other meteorological ef f ects , or 
variations in personnel effectiveness and propensity for errors. 
Random environment models conveniently postulate that environ- 
mental changes induce simple discrete-state Markovian behavior 
on the basic or primary process parameters; of most interest are para- 
meter changes that occur more slowly than do state changes in the 
basic process. 

Markov modelling of real situations usually involves simplifi- 
cations at certain crucial states. Even then, the answers to 
interesting questions may require extensive computing or simu- 
lation. Astute choices of sub-models or component models, 
e.g. the use of "phase-type" distributions for representing 
arrival and departure processes in queues can be of help, as 

23 



can the recognition (or plausible imposition) and exploitation 
of special structure; see Neuts [1981], 



24 



(6 ) Some Markov Process Problems and Models 

Here are some illustrative situations and corresponding 
Markov chain models. 

Situation (Queueing in discrete time). A servicing facility, e.g. 
a computer system or a programming (or other) consultant, or 
a communication channel, experiences single customer arrivals 
in a random fashion; arrivals enter at the discrete times 
0,1,2,3, ... only, and service completions occur only at such 
times. Discuss the nature of the delays and backlogs that 
occur. 

Model 1 . Let the probability of a single arrival at time (epoch) 
t be a. (t) , where i refers to the number present at that 
epoch. Each arrival must wait at least one time period before 
discharge, even if it immediately enters service upon arrival. 
Let d. (t) be the probability that an arrival that has been 
in service at t actually departs at t+1. Now let X (t) denote 
the number of arrivals in the system who have not yet completed 
service at time t. Model (X(t)} as a Markov chain with the 
following one-step transition probabilities: 

p. . (t) = P{X(t+l) = i + l|x(t) = i} = [1-d. (t)]a. (t) 

1 / 1 ' X i J- 

, i > l 
p. i _ 1 (t) = P{X(t+l) = i-l|x(t) = i} = d i (t) [l-a i (t) ] 

(6.1) 
p 0fl (t) = a Q (t) 



P ,0 (t) = ^V^ 



25 



P i± (t) = 1 - {[l-d i (t)]a i (t) + d i (t)[l-a i (t) ]} 

p. . (t) = otherwise; 
If the number in the system is I , so the state space is finite, 

and P i+i ft) = ° and P i i-l (t) = d i (t) theD the Probability dis- 
tribution of X(t) for any t can be obtained by numerically 
multiplying the one-step transition matrices, p_(t) , with elements 
aiven above : 

P..(t) = P{X(t) = j|x(0) = i} = element 

in ith row, jth column of 

t 

p(t) = n p(t') (6.2) 

t'=0 
This can be done especially easily in APL if the process is 

time-homogeneous, i.e. a. (t) = a. , d. (t) = d. independent of 
elapsed time. Explicit analytical solutions can rarely be 
found for non-time-homogeneous cases, let alone for time homo- 
geneous cases. If they were availabler the solutions would 
generally be very complicated and difficult to interpret. 



Model 1-A . Specialize the above to let a. (t) = a>0 and 

d. (t) = d>0. If the maximum number in the system is I, there 

is a stationary solution; put s = — , a = 1-a, d = 1-d: 

ad 

(d-a)d 



TT„ = 



dd-aas 



(d-a)s 3 (6.3) 



7T . 

dd-aas 



^ _ (d-a)as 
dd-aas 



26 



Furthermore, if s<l then the process tends to drift towards 
and even if there is no upper bound on system states the process 
is irreducible and ergodic so the long-run distribution is 

tt. = ^ s j , j = 1,2, ... (6.4) 

J dd 

a modified geometric distribution, some form of which so often 

appears in queueing problems. The above simple special case is 

well-known, but can be useful for checking the accuracy of 

computer programs used to compute numerical solutions to the 

time-dependent case • 

Model 2. Let a. h (t) be the probability that at time t there 
occurs a bunch of arrivals of size b(b = 0,1,2, ...), given that 
i are awaiting service and will not make further demands. For 
example, suppose there are I total customers, e.g. computer 
terminals accessing a central facility, and that each applies 
for service independently with probability a(t), provided that 
it is not undergoing service. Then 

a ±fb (t) = ( I ^ 1 )[a(t)] b [l-a(t)] I " i " b , 
and for k >_ , 

Pi, i+ k (t » - 5 i (t)a i,k (t) + d i (t » a i, k+ i (t) - ,6 ' 5) 



while 



Pi,i-l (t) ■ a i,0 (t,d i (t) 



27 



Again the probability distribution of X (t) can be numerically 
computed. 

Model 2 -A . Suppose the arrivals are caused by a common event 
(a "common-cause" in engineering parlance) . This might be the 
occurrence of an earthquake of large magnitude, or other en- 
vironmental shock. Let such an event occur at time t with 
probability c(t); let the probability that b arrivals 
(demands for service) occur as a result be conditionally binomially 
distributed with parameter 6. Then 

I - i \fi b n_ fi 1 I - i - b (6.6) 



^i,b (t) = cttM^e-Ei-ft] 



This can again be used to form one-step transition probabilities, 
and to calculate state probabilities at any time. The present 
model allows for a catastrophic shock situation: if 6=1 then 
all outstanding customers simultaneously demand service, i.e. 
I-i arrivals occur simultaneously. This differs from Model 2. 
Situation (Queueing with breakdown of service or preemptive 
priorities). Suppose a single server, e.g. computer facility, or 
data transmission channel, is confronted by a random arrival 
stream of basic service demands. These demands may be characterized 
by their service times, or work request durations such as the 
times required to transmit single bodies of data or digitized 
messages. In addition, these services may be effectively prolonged 
by the occurrence of interruptions, e.g. from internal server 
breakdowns resulting in temporary processor unavailability, or 



28 



from environmental noise or even intentional jamming. How is 
the queue size and waiting time of demands affected by such 
interruptions? What steps can be taken to reduce the interruption 
effects? 

Model . Nearly all classical queueing theory is most conveniently 
developed if the service times of the individual demands are 
independent"! v and identically distributed (I ID ); see however 
Jacobs [1978, 1980] for discussion of a model involving correlation 

effects. If service times are to be interrupted and repeated, 
or alternatively resumed , an interruption process that preserves 
the IID character of the basic service times allows nearly direct 
adaptation of conventional theory; such a process is one that 
requires exponentially distributed ( "memory less" ) periods 
between successive interruptions, and this will be assumed. 
Checks of the sensitivity of results to this reasonable 
assumption can be made by simulation. 

If interruptions of IID duration X (df F (x) ) occur at 
IID Expon(A ) intervals, then the time to complete the ith 
basic service (low-priority) is, provided service can resume 
after each service 



C i = S i + X l + X 2 + •'• + X I(S.) ' (6 ' 7) 

where S. is the ith basic or low priority service time 
1 — 

(df F (x) ) , X. is the duration of the jth interruption, and 
I(S.) is the number of interruptions that occur during S.. 



29 



Given S., I(S.) is Poisson (A ) , and the Laplace-Stielt jes 

11 n 

transform of C. is, in terms of 



F c (s) = F s [s+A H (1-F x (s) }] 



(6.8) 



and hence 



E[C] = E[S](1+A H E[X]} 



E[C 2 ] = E[S 2 ]{1+A U E[X]} 2 + E[S]A U E[X 2 ]. 

ri n 



(6.9) 



If, on the other hand/ basic services that are interrupted must 
begin again from scratch, i.e. services must repeat , then the 
i th completion time becomes 

C. = S. + X 1+ S| + X 2 + S* + ... + X I(S _ } + SJ. (Si , (6.10 

where S '. is the ith interrupted basic service time that must 
D — 

be repeated. The L.-S. transform is 



F (s) = E< 
c 



-(A H +s)S 



A H +s 



A TT r - (s + A rT )S 



1-e 



H 



■>, 



6.11 



and by differentiation of the latter expression, 

r AS- 



E[C] = ^E 



H 



1 E[X] + \ 



II 



(6.12 



t A..S 2-, 



E[C 2 ] = 2E (e H -l) 



E[X] + i 

J I A 



11* 



+ 2E 



Se 



(E[X] + \- 



A 



H 



+ 



V 



1 2 



E(e " )-l) (E[X"] + 2E[x].f- + ^ 



30 



In order to assess queueinc delay, look at the process, {N,,d=0,l 2 
describing the number of basic demands at the server at the instants 
just following departures; { N d , d - 1,2, ...} is an embedded Markov 
chain provided basic arrivals are (compound) Poisson. The 



N,. 1 = H , + A(C,.,) if N , = 
d+1 d d+1 d 



= N + A(C,, 1 )-l if N, > 1 . 
d d+1 d — 



(6.13) 



Here H is the number of basic (low-priority) demands made at 

the beginning of a basic service busy period initiated by the 

appearance of a high-priority demand, and A(C,) is the number 

of basic demands made during the dth basic completion time. 

Express an arbitrary H as follows 

A L 
with probability t — t-t — 

H - L H 

^ H (6.14) 

A(X) with probability t — -^ — 

A L +A H 

It follows by conditional expectations that the embedded chain 

is ergodic if E[A(C)] < 1, and that then the long-run probability 

of system emptiness at an embedded time point is 

1 - E[A(C)3 



p = lim P{N =0} = - r ^ J (6.15) 

O j d ELHJ 

and the long-run expected occupancy is 

E[N] = 2(1-E[A(C)J) { E[(H + A ( C )) 2 ]P + E[(A(C)-1) 2 ](1-P )}. 

(6, 16) 
Delay can then be estimated by use of Little's formula. A version 

of the above formulas correct in continuous time may be found using 



11 



results of Gaver [1962],; the difference between embedded and con- 
tinuous time becomes comparatively negligible if the basic 
traffic intensity E[A(C)] is close to, but below unity ("heavy 
traffic"). Since computer system monitoring devices sample 
system state at the moment an event occurs (e.g. at a departure 
instant) a theoretical account of the queue at such moments 
(imbedded times) is of direct interest. 

An alternative approach to the long-run distribution of 
delay of an arriving basic demand is by way of Wald's identity 
or martingales; see Feller [1966] . This will actually handle 
waiting-times when the basic service inter-arrival intervals are 
IID, but otherwise arbitrarily distributed. Still another ap- 
proach is via the Takacs-type integro-dif f erential equation^ see 
Kleinrock [1976] for an account. 



32 



The previous situations have been discussed in terms of long- 
run probabilities. Frequently questions involving the time until 
system failure (or restitution to operational condition) are more 
important. 

Situation (Redundant repairable systems) . A particular system 
function can be performed if at least k out of n system com- 
ponents ("machines") function. For example, electric power is 
available if at least one generator is working out of two that 
are installed. Suppose that the system components are all operative 
initially but fail randomly; failures are immediately detected, 
but repairs are of random duration, so several machines can be 
down, all awaiting repair completion. Question: how long will 
it be until I = n-k+1 components are simultaneously in a failed 
condition? The time until this occurrence is the time to failure 
of a k - out-of - n system. 

Model 1 . It is now convenient (but possibly unrealistic!) to assume 
that the machines fail independently and after exponentially dis- 
tributed times in operation, each with rate A. This simplification 
may be relaxed, but at the price of expanding the state space. 
Assume too that the repairs occur in a Markovian manner, e.g. (but 
not necessarily) at rate y -Min (N (t) ,R (t) ) , where N(t) is the 
number of machines failed and down for repair at t, and R(t) is the 
number of repairmen on duty. This is the classical machine repairman 
problem; see Feller [1966], and Cox and Smith [1962]. Usually R(t)=r, 
a constant, although provisions may be made for automatically in- 
creasing repair effort when redundancy reserves become dangerously 
low. In other words, N(t) is a simple birth-and-death Markov 
process, wherein jumps in state, N(t), occur at exponentially dis- 
tributed intervals or sojourn times , S^ for the sojourn in state i, 

33 



transitioning always to neighboring states. In the present 
situation as A -* 

P{N(t+A) = i + l|N(t) = i} = A(n-i)A + o(A) 

= A. A + o(A) 
1 

P{N(t+A) = i-l|N(t) = i} = y min(i,r)A + o(A) 



(6.17) 



= y . A + o (A), 
abbreviating the general transition rates, to A. and p. . 

In order for system state to reach i+1 from i for the 

first time it must either do so on the first transition out of 
state i, or else drop back to i-1, return to i and try again 
Thus if U. is the local first passage time from i to i+1: 



U. = inf{t: N(t) = i+l|N(t) = i}. 



(6.18) 



write 



U. = S. + { 
l l 



with probability p. . , 



U! , + U! with probability p. . , 
l-l i r J r i f i-l 



A .+y . 
l l 



(6.19) 
y i 



A .+y . 
l l 



where U! has the same distribution as U-. The above repre- 
sentation allows immediate derivation of the Laplace-Stielt jes 
transform of U. by conditional expectations. The result is 



-sU 



A . 

E ^i (s) = s+A.+y.[l-ij,. I (s)J ' 1 = 1 
l i L ^1-1 



= 12 

(6.20) 



^o (s) = A-Ti • 



34 



Furthermore, since the first-passage time from N(t) = i to j > i 
can, on the basis of Markovian assumptions, be expressed as 



T. . =U. + U. , , + 
i] 1 i+l 



+ U 



j-l 



(6.21) 



where {U. , , k = 0,1, ...} are independent rv, the L.-S 

transform of T. . is 

ID 



-sT. . 



j-l 

n 

k=i 



n ij> k (s) 



(6.22) 



and the cumulants (moment-like quantities; see Cramer [1946]) 
can be expressed in terms of those of the U- . Here are a 
few moments of the U.; recursively expressed and hence easily 



computed : 



:[u ± ] = b{ 1 + ^i E[u i-i ] }' 



K, + UiEEu.^]} + Jr E t u i-i3 



AT 

l 



E[u 3] . 6 (, + p.ECu.^]} 3 + ^ jl + M i E[0.. 1 ]}E[u2_ 1 ] 



and 



+ £ K[UJ] 

1 



(6.23) 



E[nJ] - % {i + v ± EtUi.^} 4 + ^ {l + ^[U.^]}^^] 

A . A . 



6y 



2 8y . 



E ^ u i-i ] + -TT i 1 + ^i E f u i-i ] E[u i-i ] 



i „r„4 



+ ji eCU^] 



35 



From these, standard variance, skewness , and kurtosis measures 

can be easily computed. For the repairman model discussed 

initially it can be shown that if the expected time to system 

failure, E[T n ] , is "long" then T /E[T „] resembles an 
L o£ ox/ L ol 

Expon(l) rv. 

Model 2 . (Catastrophic failures) . Suppose that in addition to 
the independent random failures there is a catastrophic event 
that "kills all operative machines simultaneously; let it occur 
after a time C ^ Expon(v). Then the system failure time T* 
has distribution given by 

P^T* £ > tl = p{t o£ > tl e" Vt (6.2.) 



and from this 



E[e ST °^ = iz£ 



- (s+v) T 



s+v 



(6.25) 



from which moments can be generated; see Chu and Gaver [1977], 
It is sobering to note that if v , the mean time to catastrophe 
occurrence, becomes small or even comparable to e[t „ ] , then 
the mean time to redundant system failure is essentially v , 
and redundancy alone may not improve system reliability. 

Model 3 (Simultaneous repair) . If the system is not under con- 
stant surveillance, but instead is inspected at random times 
(rate y) and then repaired in negligible time, the number of 
down machines at t may jump essentially instantaneously, either 
to zero (perfect and rapid repair) , or to some lower point 
(imperfect repair) . In this case the basic "skip-free up" 



36 



character is retained, but now 

with probability p 

U. = S. + < 
1 l 



i,i+l 



U! + U! + U! + ... + U! . with probability p. 

j = 0,1,2, . . . i. (6.26) 

Note that if j = then repair is completely ineffective 
Conditional expectations now give the L.-S. transform 



-sU. 
E[e 1 ] = ^.(s) = 



V S)P i,i+l 



(6.27) 



l-n ± (s) I 



n i> 

=l l r=l 



l-r 



i/i-: 



-sS. 



where here n. (s) = E[e x ] = a. (a.+s) . To specialize this to a 



l i 



repair model, introduce a Binomial distribution for successful repairs 



a . = X (n-i) + u 



p i,i+l = X(n " i) a ^ 



-1 

i 



(6.28) 



p . . . = pa . 



l 



P j (1-p) 



i-D 



j = 0,1,2, . . . , i 



where p represents the probability that an individual down 
machine is indeed repaired just after inspection (neglect the 
duration of repair times) . The Binomial model assumes that 
repair success is independent across machines, which may be 
inappropriate in case similar causes give rise to the failures. 
Differentiation or direct expectations yield moments of U. 
and eventually of T p . 



37 



(7) Diffusion and Fluid Approximation 

While classical discrete state space Markov process ideas 
can often be used to model some quite interesting situations, 
the analytical results obtained frequently emerge only in terms 
of incomprehensible transforms, or in other somewhat obscure 
form. Not infrequently the difficulty that induces complexity 
can be traced back to the influence of boundaries upon the 
process transitions. If one examines a rather heavily loaded 
or congested service system, however, it is apparent, first, 
that the state changes may appear almost negligibly small relative 
to the system state magnitude (e.g. length of queue) itself, 
and, second, that annoying boundaries, particularly that at 

zero, are visited infrequently although their influence 

may still be crucial. These remarks hold true not only for 
simple one-dimensional processes, such as those used to describe 
congestion at a single servicing facility, but also for much 
more complex situations involving the interaction of several 
servicing processes. 

An attractive approach to problems involving many customer 
arrivals occurring rapidly and generating considerable queueing 
is, then, to treat them by the method of diffusion process 
approximation. For details concerning the rigorous details of 
diffusion mathematics see Feller [1966]; in brief summary recall 
that a diffusion is a possibly vector-valued Markov process on 
the real numbers that typically moves continuously, governed by 
a drift (infinitesimal mean) and a diffusion (infinitesimal 
variance) parameter. 



38 



This approximation has been employed by Kobayashi [ ] for certain 
cyclic networks of queues such as are encountered in multipro- 
gramming computer systems. See also Gaver and Shedler [1973], 
and the important work of M. Reiman [1982], and particularly 
G. Newell [1979] and also McNeil and Schach [ ] . In this section 
the use of diffusion will be briefly illustrated, and some ex- 
perience with the results will be recounted. 

Situation (Waiting time or backlog at one server) . Suppose a 
single servicing facility is confronted by random arrivals that 
bring with them contributions to work load, expressed as re- 
quired processing times. If the facility processes them in 
order of arrival, what is the backlog at time t? The facility 
is heavily loaded, so that it is seldom idle. It never turns 
away customers, i.e. infinite buffering is possible, nor do long 
delays discourage those waiting, causing defections or balking. 

Model . Assume that arrivals occur in a Poisson (A) process, 

and that the generic processing time S has df G (y) ; successive 

processing times are IID. This model has been studied by Takacs 

] who derived an integro^dif ferential equation for the df 
of backlog or virtual waiting time W(t). Although the formal 
solution of that equation can be obtained, it is in a somewhat 
complicated form, not conducive to immediate insights. It is 
tempting to take an alternative, somewhat heuristic approach. 
Intuitively, if p = AE[s] (= expected total load increment per 
unit time) > 1 (= processing or output rate per unit time) , the 
backlog grows at rate p - 1 > 0. Furthermore, the backlog process, 



39 



W(t), "ignores the boundary" at W = after a time, and 
eventually w(t) appears approximately Normal/Gaussian over an 
interval (t, t+A) , with mean (p-l)A, and variance Ae[s 2 ]A. 
Importantly, also, the process seems to grow by accumulating in- 
dependent, nearly Gaussian, increments. 

If p<l some difficulties occur because W=0 is an im- 
permeable reflecting boundary, but the basic scenario is still 
the same: if p is close to unity W(t) moves in nearly Gaussian 
increments, but occasionally interacts with the boundary at 
zero. Question: what is the long-run behavior of the delay 
in such a process? 

The Takacs (or forward Kolmogorov) equation for the "exact" 
process is 

- H + H = -AF + A / F(x-y,t)G (dy) (7.1) 

o 

where F(x,t) is the df of W(t); initial and boundary conditions 

are necessary but are suppressed. 

If F(x,t) is only appreciable when x is large, and if the 
magnitude of a typical S is also small compared to x then 
it becomes plausible to Taylor-expand the F (x-y , t) -term to 
three terms and integrate; the result is 

|| = (1 . AE[s]) | + m|!ii!| 

which is the well-known forward partial differential equation for 
Brownian motion with drift. Impose the boundary condition that 
F(x,t) = for x < 0, and the equation for F, an approximation 



40 



to F, the "true" distribution of w (t) r emerges ; the latter 
equation can be explicitly solved in terms of error functions 
and exponentials (see Newell[l97l]) for all t; i.e. the transient 
solution is actually readily available. With some further 
effort one can impose an upper boundary at x > to repre- 
sent a finite buffer size. The approximate steady-state 
(t -*■ °°) distribution turns out to be 



F (x) = exp 



-2(1-P) x 
AE[S 2 ] 



p < 1 (7.3) 



which often is, for large p, usefully close to the behavior of 
F (x) , the long-run solution to the Takacs equation. 

Note that the parameters of the above differential equation 
can be obtained by equatina infinitesimal mean and variance of the 
assumed Poisson arrival process to the corresponding quantities 
for the diffusion. It is interesting that, when available, a 
martingale approach to the problem, of Gaver and Shedler [1973], 
produces a different exponent that yields better approximations 
for moderate traffic intensities, especially if the service time 
distribution is very long tailed (more skewed than the expon- 
ential) . 

Turn now to a more complex example, involving the interaction 
of two traffic dreams. 

Situation. At a node of a communication network there are a 
total of c+v channels (servers) ; voice messages are exclusively 
assigned to the v channels, and data messages are assigned 
to the c channels, but may also utilize any unused voice 



41 



channel capacity under the stipulation that voice has preemptive 
priority and may displace any data encroaching on its (v-channel) 
territory. Question: what is the nature of the delay experi- 
enced by the data, which is allowed to queue up in a buffer in- 
definitely? 

Model (Markov assumptions). Voice traffic is Poisson (A) with 
Expon (y) service times; voice is a loss system, so immediately 
the steady-state voice loss rate can be calculated using the 
Erlang-B formula. The data, however, operates in a 
random service environment modified by 

voice needs. Data arrives in an independent Poisson (5) process, 
with independent Expon (n) service times. Typically 6 >> A, 
and n >> y. Data, with state variable X (t) , is an M/M/S 
system where S = c+v - V(t), V(t) being the number of voice 
messages in service. Clearly (X(t), V(t)} is a bivariate 
Markov process, but one difficult to analyze exactly; see 
references in Lehoczky and Gaver [1981] for other approaches 
to the analysis. 

Now typically n/y (.= Data arrival rate per voice service 

4 
time) is very large, possibly 10 . Furthermore often 

p, 5 6/ri > c, so some voice channel usage by data is necessary 

in order that all data be handled and there is not an evergrowing 

queue. The appropriate traffic intensity parameter for the 

system is seen to be 

P = [p d + P v (l-q)](c+v) _1 (7.4) 



4? 



P V /v! 
v 



where p = A/y, and q = — 

v v . 

I pVj: 

j-0 v 

is the probability that the voice system rejects an arriving 
voice message. If p becomes large under such circumstances, 
Gaver and Lehoczky [1982] show that X (t) behaves like a Wiener 
process with reflecting boundary, precisely as was mentioned in 
the previous example. Actually Lehoczky and Gaver [19 82] 
assumes that data input acts as a fluid , with no variability. 

By further, more intricate, methods involving the convergence 
of semigroups of operators developed by Burman [1979], it is 
shown in Lehoczky and Gaver [19 81] that the long-run distribution 
of X (t) is exponential with mean 



where 



v-1 



p J + — 
d yp 



- I (Tt/TT ) 

r i=0 
(c+v) (1-p) 



(7.5) 



and 



tt. = pV(i!)/ 
1 v 



r v 



I P^/(iD 



L i=0 



T k = I 7T i (i-p v (l-q)) 



i=0 

Numerical work indicates that the diffusion approximation is 

reasonably accurate if the traffic intensity is quite high, say 

if p > 0.95; otherwise, for smaller s, the accuracy is not as high 

It would not be surprising if a refined method for fitting drift 

and diffusion coefficients would lead to improved results. 

Finally, the difference between the refined treatment of 

basic data as a Poisson, and the simpler treatment by a 

fluid (yielding a model that can be solved exactly) resides 

43 



merely in the addition of the term p , in the numerator of the 
expression for the mean. 

The diffusion approximation can be utilized to evaluate 
another interesting measure of system effectiveness, namely 
the expected total waiting time, in job or data-packet hours, 
expended during a busy period for data traffic that starts with 
x present. 

Model . Let A(x) be the expected total waiting time during 
a busy period when the initial number of jobs is X(0) = x > 0. 
Condition on process change, Z (A) , during the initial short time 
period (0,A) to get 

e|/ X(f )dt' |X (0) = x, X(A) = x + Z (A) | E 



(7.6) 



A(x;Z (A)) = xA + A(x+Z (A)) + o (A) . 
Now Taylor-expand and remove the condition on Z(A): 



2 

A(x) = xA + A(x) + A (x)y(x)A + A (x) • °—^- A + O (A) (7.7) 

or, collecting terms in A and letting A ■*■ , 

= x + y(x)A x (x) + |a^ x) A xx (x), (7.8) 

to be solved subject to A(0) = 0- clearly restrictions on y (x) 

are necessary in order that A(x) be finite. 

2 ?. 

I* y (x) = y < a ' (x) = a , both independent of x, the eauation 

can be solved directly to give 



AA 



A(x) =h* 2 + vi * ' <7 - q) 

y<0 when traffic intensity p<l. A similar expression 
for compound Poisson (A) inputs is 

A(x) = r7 f- T + X E[A2] , p<l, (7.10) 

' u P; 2(l-p) 

where A represents the number of items (packets) arising at 
a single request; p = Ae[a]. The moral is that variability 
(measured by a or E[A ]) can greatly increase expected 
total waiting time, particularly when system loading is high 
(p close to unity) . 

For application of this "area under a random path" to discussing 
total wait during a road traffic jam see Gaver [1969J . See also 
McNeil [1970] for generalizations. The same backward argument is 

well-adapted also to studying problems of optimum investment 
decisions. 



45 



8 . Renewal-Theoretic Modeling 

Ideas of renewal theory and recurrent events are extremely 
useful for many purposes in stochastic modeling. Recognition 
of the occurrence of one or more recurrent events or "renewal 
points" in the development of a model process points the way to 
writing down simple forward or backward type equations for 
probabilities or expectations. Frequently analytical information 
can be extracted from such equations, particularly that relevant 
to long-time or other asymptotic process behavior. If more in- 
formation is desired it can be obtained by use of transform 
techniques, by numerical computation, or by Monte Carlo simulation 

Mathematical definitions and properties of renewal processes 
are well presented by Feller [1966], Karlin and Taylor [1975], 

and Cox [1962], among others. There follow a few situations and 
suggested models based on renewal theory that illustrate the 
basic notions. We also comment on the relevance of the models 
and results obtained to real situations. 

Situation . A machine, e.g. computer system or component thereof, 
or human operator, etc. , operates properly for a period of time, 
fails, is restored (or restores itself) to service and operates 
properly again for a different time, fails again, and so on. 
Questions: how many failures are likely to occur in a given 
fixed period of time, say a year? The answer to such a question 
will help to guide decisions concerning logistics (necessary 
spare parts) and employment of repair personnel. How long 
a time will elapse until the k — failure? Suppose the times 
to restore failures vary; what is the likelihood that a "chance" 



46 



user of the system will find the machine down for repair when 
he or she needs it, and how long will the wait to service 
restoration last? These are only a few of the many questions 
that might be asked. 

Model . The classical renewal theory model for the situation 
described portulates that times between successive failures, 
{X-, i = 1,2, ...} are IID positive rv with distribution 

F (x) . That is X, is the time until the first event (here 

• -, st , • th 
failure) , and X. is the elapsed time between the l-l — and 1 — 

event. X. can be either a discrete random variable, so failures 
occur at regular intervals, say hourly, a continuous random 
variable, or a mixture. For the moment assume repairs to take 
a negligible time. 

Mathematical results for this model are simplest and 
nicest when the IID assumption is fully exploited. Under that 
assumption (and even more generally) the counting process, N(t), 
giving the number of renewal events (failures in time t) has 
probability distribution 

P{N(t) = n) = f£* (t) - F^ n+1) *<t) (8.1) 

where * refers to convolution. For long time (t -*- °°) and under 
suitable mathematical restrictions 

M(t) = E[N(t)] ~ ^J (8 ' 2) 



Var[N(t>] - Var[X] , t 
" (E[X]) 3 



47 



and N(t) ^ Normal, with the above parameters. Of course exact 
analytical solutions can be obtained if sympathetic distributional 
models are assumed: taking X ^ Expon. yields the Poisson 
distribution, and X ^ Gamma or Erlang also produces rather neat 
closed-form solutions. Any discrete-time distribution for X 
can be numerically convolved, conveniently using APL. This 
helps to answer questions about the number of events in (0, t) , 
provided the IID assumption is palatable. If finite-time results 
are needed resort can be made to numerical summation, using a dis- 
crete time model, to approximation by a standard, tractable, 
distribution such as the Gamma followed by transform inversion, 
or by simulation. 

In what follows we illustrate the diverse utility of 
backward conditioning arguments, leading to renewal integral 
equations. 



48 



Situation . (Incorrect repair possibly due to human error) . Each 
time a repair is made there is the chance that it will be incorrect, 
and that the subsequent time to failure will be short. Suppose 
that incorrect repairs tend to bunch together in runs; describe 
the number of failures that occur over a (long) time period. 

Notice that a similar situation describes a clustering 
scheme of arrivals to a repair facility or a communications center. 
Model . Assume that the generic time to failure of the system is 
X when repair is made properly, and X' when repair is incorrect; 
the intermingled sequence of X and X 1 quantities are condi- 
tionally independent. Furthermore, let the probability of a 
correct repair at failure n be a, 5 a 5 1, if the repair 
at the time of previous failure n - 1 was correct, and 3 = 1-6 
< 3 5 1 (a^3) if it was incorrect; the sequence of correct 
and incorrect repairs is thus modeled as a stationary 
ergodic Markov chain. This, of course, does not represent system- 
atic improvements in repair capability, although a transient 
chain could serve for that purpose. 

Let M(t) (M'(t)) denote the mean or expected number of 
repairs in t, given that the first repair was correct (incorrect) ; 
think of the first repair (manufacture) as occurring at t = 0. 
Argue that 

1 if X > t; 
M(t) = \l + M(t-X) if X<t and the 2nd repair is correct; (8.3) 
1 + M' (t-X) if X<t and the 2nd repair is incorrect. 



Likewise , 



49 



1 if X' > t; 
M' (t) ]1 + M(t-X') if X' < t and the 2nd repair is correct; (8.4) 
.1 + M* (t-X 1 ) if X' < t and the 2nd repair is incorrect. 



Now if the various conditions are removed then according to the 

model there results the two linked convolution integral equations 

t t 



M(t) = 1 + a 



M(t-x)F x (dx) + (1-a) M(t-x)F x , (dx) . 



and 



(P. 5) 



M' (t) =1 + 3 



M(t-x)F x , (dx) + (1-3) 



M(t-x)F x , (dx) 



(8.6) 



In turn, these two equations are susceptible to transforming: 



multiply by e 



-st 



and integrate to get 



M(s) = - + aM(s)F v (s) + ctM 1 (s)F v (s) 
S X X 



M (s) = 



± + 3M(s)F x ,(s) + 3M(s)F x ,(s); 



matrix notation is natural here, especially if more than two 
repair states are used. If one then solves and collects terms 
to order (1/s) as s -> , Tauberian theorems, cf. Feller [1966] 
show that for large t 



M(t) - M 1 (t) ~ 



a + R 



aE[X' ] + 3E[X] 



(8.7) 



and a little reflection shows that this is entirely sensible. 

In similar ways variances can be written down, and an approximate 

Normal/Gaussian distribution for total failures may be derived. 



50 



The model can also be extended to account for the existence of 
non-zero repair times, and total availability studied. 

Here is another, somewhat more complicated, application 
of renewal theory ideas, now to an inspected system problem. 
Situation . (Standby system availability) . A system, such as an 
emergency electric power source, is usually in a quiescent or 
cold standby status, but occasionally is called upon to fill a 
function (e.g. generate power) . Various systematic plans might 
be devised for assuring reasonably high system availability, or 
probability of satisfying a demand when one occurs. One such 
plan is to inspect infrequently so long as no failure is detected 
between inspections, and otherwise inspect more frequently until 
evidence of need seems gone. Problem: develop a model to evalu- 
ate such an inspection scheme. 

Model . Let the inspection plan be to inspect at long intervals, 
{L., i=l,2,...} until such time as an inspection reveals a fail- 
ure, and then switch to short intervals {S., i=l,2,...}, con- 
tinuing until there has been a run of r(10,say) failure-free 
short-interval inspections, at which time switch back to long 
intervals; continue indefinitely. A measure of effectiveness is 
the long-run point availability of the system, i.e. the proba- 
bility that the system is failure-free on the occasion of a 
demand. 

To evaluate such a rule, allow the L. and S. sequences 
to be IID and mutually independent, with dfs f t/ x ) and F S ^ x ^ 
respectively; if desired these latter can be specialized to con- 
centrate at fixed values (e.g. 14 days and 1 day, respectively. 



51 



Furthermore, let A be the failure rate of a system failing at 
exponentially distributed intervals even when "cold." In order 
to analyze the system by renewal theory it is worthwhile to look 
at the periods during which inspection are infrequent, called 
L-eras , and those alternating with them, during which inspec- 
tions occur frequently, called S-eras. Note that a demand can 
occur during either type of era, L- , and S-eras constitute 
an alternating renewal process. 

The analysis involves the following components. 
(A) Distribution (density) of L-era duration. 

Suppose an inspection has just been completed at t = 
and nothing amiss has been detected. Furthermore, suppose that 
this inspection marked the end of the previous S-era, so an L-era 
is just beginning. Let a (dt) be the probability that the 
present L-era will last for time (t,t+dt), or, loosely, until 
exactly t. One can now write down a renewal equation for 



a L Cdt) : 



a L Cdt) = (l-e" Xt )F L (dt) + 



t 
e" At 'F T (dt')a_ (dt-f) ; (8.8) 



the first term on the rhs means that the L-era terminates with 
the first inspection, meaning that the unit has failed before 
L, . The second term represents survival through the first inspec- 
tion at which time the process renews itself or starts anew; final 
failure occurs at time t-t' thereafter. Failure or no failure 
at first inspection are mutually exclusive and exhaustive events, 
and so the result is a renewal equation for a (dt) . Introduce 
transforms to find 



52 



F (s) -F (A+s) 

a (s) = E[e -] = -±! £ . (8.9) 

1 -F L (A + s) 



The mean of the L-era duration, denoted by L, is 



E [L] = E[L] . (8.10) 

1-F L (A) 



(B) Distribution of S-era duration. 

If an inspection has just revealed a failure, an S-era 
begins immediately (take inspection and repair to be instantane- 
ous for the moment). Let a q (dt) be the probability that an 
S-era lasts for time t. Then the following renewal equation 
may be written: 

t 
a g (dt) =e" At F^*(dt) + h(dt')a s (dt-t'); (8.11) 



the auxiliary function h represents the probability that an 

inspection reveals a failure before the termination of the S-era 

in progress; this causes the frequent inspection to start over, 

i.e. starts the S-era afresh. In terms of transforms, S denoting 

an S-era duration, 

_ (F c (s+X)) r 

a c (s) b E[e" S -] = — 2 , (R.12) 

b l-h(s) 



and 



53 



h(s) = [F s (s)-F s (X + s) ] •{! +F s (s + X) + (F s (s+X)) 2 + ... + (F s (s+X)) r 1 } 



(8.13) 



F (s) -F (s+X) 

-S 2 [1-(F fs+X)) r ] ; 

1-F (s+X) b 



The latter transform expresses the probability that the necessary 
run of r successes is interrupted by a failure, and must begin 
again. From the transform comes the expected length of an S-era: 

(F_(X))" r -l 

E[S] = E[S]{ — } . (8.14) 

1-F S (X) 



(C) Probability that system is available during an L-era. 

The probability A (t) that the system is up at time t 
after the beginning of an L-era is simply e , the probability 
of no failure, and the transform is 



e St A L (t)dt = A L (s) = (X+s)" 1 . (8.15) 



(D) Probability that system is available during an S-era. 

If A (t) is the probability that the system is up after 
an S-era has progressed for time t, by renewal 

t 



A s (t) = e At F s r *(t) + 



h(dt')A (t-f) (8.16) 







where h has appeared before, under (B) . It follows that 



54 



A s (s) = 



1 - (F s (s+X)) r 
(s+X) [l-h(s) ] 



(8.17) 



(E) Overall availability at t. 

Let A(t) denote the availability of the system at t. 
Again by backward renewal argument, and starting at the beginning 
of an L-era, 



A(t) = A L (t) + 



A s (t-f)a L (df) + 



A(t-t')a L * a s (df ) ; (R.18) 



transforms then give 



A(s) = 



A L (s) +A s (s)a L (s) 
1 - a L (s)a s (s) 



(8.19) 



a Tauberian theorem now shows that 



lim A(t) = 



A L (0) +A S (0) 
E[L] + E[S] 



so long-run point availability is 



A(oo) = (A 1 ) { 



(F S (X)) 



-r 



E[L] (1-F L (X)) 1 + E[S] ((F S (X)) r -l) (1-F S (X)) 1 



} . (8.20) 



The expression for A(°°) is easily evaluated numerically 
if expressions for L-interval and S-interval transforms are 
obtainable. It is sometimes useful to interpret 



55 



A (s) • s = 



A(t)e St sdt (8.2,1) 



as the availability of the system upon demand , the demand now 
occurring at a random time D, D having the exponential 
distribution with mean s . In this case the initial condi- 
tions matter (they do not, in the long run) , and a different 
availability figure is obtained depending upon whether the 
system is initially in an L-era or an S-era. 



56 



4. Additional Modeling Topics 

In this section are outlines of certain modeling topics 
that formed the basis for cooperative research at ULB . 



57 



(1) Distributional Sculpturing or Inverse Modification 

Standard distributions such as the Exponential or 
Gamma in particular, and also the Normal, log-Normal, and many 
others may reasonably and conveniently serve as components of 
a stochastic model or be used to summarize data distributions. 
For instance, inter-arrival times to service systems may appear 
approximately Exponential, service times nearly Gamma or log- 
Normal, and so on. On the other hand, a systematic departure 
from such a standard may be revealed by an analyzing actual 
data. It is then frequently possible to alter conventional 
distributions, or, equivalently , transform the random variable 
in simple ways in order to represent empirical reality more 
closely. Here are two conventional examples; there then follow 
some more general procedures. 

Example 1 . The Weibull distribution is often utilized to 
represent times to failure of components or times between sys- 
tem demands. Let T be a random variable (time, for instance) , 
then T is distributed in a Weibull manner if 



( 



F T (t;a,3) = P{T < t} = ■ 



1 - e 



-at 1 



, t>0,a>0, 3>0 
t < ; 



(1.1) 



The Weibull density is 



3 



f T (t;a,3) = e at " a 3 t 6 " 1 



(1.2) 



58 



Furthermore, the Weibull hazard or failure rate at age t is 

~ P{T 6 (dt) |T > t} E k T (t;a,3) 

f (t;a,3) R , 

= aBt P L ; (1.3) 



1 - F T (t;a,B) 



The latter formula and also (1.2) imply that if 3=1 the 
Weibull is actually an Exponential; for this case "age" or 
"time in service" as measured by t does not influence the 
probability of failure in the next short time interval 
(t,t+dt) . On the other hand 3 > 1 implies that the rate of 
failure having survived to t — the hazard--increases with age, 
while, 6 < 1 implies that hazard decreases with age. If 
3 < 1 the Weibull right tail is longer than that of the 
Exponential (extremely positively skewed) , while if 3 > 1 
the positive skewness is less pronounced. The Weibull r.v. , T, 
is actually only an Exponential r.v., X, transformed or 
disguised, for if 

P{T < t} = 1 - e' at = P{T 6 < t 6 } , (1.4) 



then 



~ r m 3 i -1 -ax 
P{T < x} = 1 - e 



Consequently, T = X , an Exponential r.v. , and T = X is 
a representation of a Weibull in terms of an Exponential. Sup- 
posing that one wishes to simulate a Weibull (a, 3) r.v., then 
one merely simulates an Expon (1) r.v. and raises it to the 



59 



(1/3) power, later multiplying by a~ . Since the power 
transformation is monotonic increasing, the quantiles of 
Weibull and Exponential are related by 



t(p) = (x(p)) e (1.5) 



Example 2 . The log-Normal distribution is a favorite model 
for system repair times (see Kline and Alvooq [l Q nn]« it has r^any other 
uses, and even some rational persuasions for its apparent re- 
semblance to empirical distributions (see Aitchison and Brown 
[1957]). Say that X is log-Normal if in x = Y is Normal 

(it doesn't matter what the base of the logs is.'); specifically 

2 
Y ~ N(\i,o ). Here are some properties: 

. ^1.2 2 

• m k = e = E[X K ] 

u+~-a _ 2 9 9 

• n^ = E[X] = e , Var[X] = (E[X]) / (e° -1) = (E[X]) n 

Y 1 (X) = skew[X] = n + 3n (1.6) 



y 2 (x) =kur[x] = n 8 + 6n 6 + I5n 4 + I6n 2 



2 

Median[X] = e y , Mode[X] = e y ° 



Supposing that one wishes to simulate a log-Normal random 

2 
variable, one simply simulates a W(y,a ) r.v. , Y, and 

Y 
exponentiates X = e 



60 



The above two familiar examples illustrate the for- 
mation of new and useful random variables and distributions 
by simple transformation. An intuitively appealing way of 
looking at certain transformations is in terms of modifica- 
tions to familiar random variables or quantiles by convenient 
shaping factors . This process will be called 

Distributional Sculpturing: Let X be a basic r.v., for 



example an Exponential, Normal or log-Normal. Define 



Y = Xs(X) 



(1.7) 



where the shaping function s(X) is designed to conveniently 
convert the basic r.v. X into a shaped version, Y , having 
desired distributional properties. Some examples of the prop- 
erties often desired in practice, along with suitable — but not 
unique--shaping functions, now follow. 



(A) Skewness-Producing Shapers 

Example: X is a positive basic r.v., e.g., Exponential. 

(i) s(X) = 1 + AX £ , A > , I > d.8) 



(ii) s(X) = e 



AX' 



Effect: 



Y = X s(X) (i) 

(ii) 



X 



if X "small" 



X + AX 

,1 



£+1 



>> X 



Xe 



AX' 



much greater 
than X if X 
"large" 



61 



The particular shaping functions (i) and (ii) both leave small 
values of X unchanged, but considerably expand large values, 
thus transforming the distribution of X , e.g., the Exponen- 
tial, to one that is nearly Exponential near the origin, but 
having a relatively long right tail. For example, take shap- 
ing function (2) with 1 = 1 and apply to X Exponential. 
The shaped distribution becomes 



F Y (y;A) = 1 - exp[- ( /1+4 ^ ^) ] . (1.9) 

Examination shows that for small y (Taylor series) the dis- 
tribution of Y is nearly unit Exponential, while for large 
Y it resembles a Weibull with shape parameter (3 = 1/2. Shap- 
ing function (ii) has an even more pronounced effect on the 
right tail. Note that both shaping functions (i) and (ii) 
yield monotonic increasing transformations from X ■* Y , and 
that given by (i) is sometimes explicitly algebraically in- 
vertible (solve quadratic equation when 1=1, cubic when 
1=2, quartic when £ = 3) , while that of (ii) is not. Also 
note that useful transformations result when parameter A < 0: 
this actually may result in right tail truncation (severe 
shortening; such transformations are no longer monotonic. 
Moments . 

The moments of shaped or sculptured r.v.s. can some- 
times be conveniently calculated. For the present represen- 
tations: 



62 



(i) Y = X(l + AX^) ; 



m (Y) = m 1 (X) + Am 1+Jl CX) 



m 2 (Y) = E[(X(1 + AX £ )) 2 ] 



= m 2 (X ) + 2Am 1+£ (X) + & 2 ™ 2 +2l {X) (1.10) 



Var[Y] = Var[X] + 2ACov[X, X £+1 ] + A 2 Var [X £+1 ] . 

Also, to indicate the dependence between the stretched r.v. 
Y and the basic r.v. X , 



m k (Y-X) = A k m 1 CX U+1)k ) = A k m u+1)k (X) 



so 



and 



Var[Y-X] = A 2 Var[X £+1 ] (1.11) 



Cov[Y,X] = Var[X] + ACov[X f X £+1 ] . (1.12) 



The quantiles are directly and simply related by 



y(p) = x(p) (1 + A(x(p)) £ ) . < p < 1 (1.13) 

I 
AX 
(ii) Y = Xe ; but for the present consider only 

£ = 1 . Then the k moment is expressible in terms of the 

derivatives of the Laplace transform of X . Note that the 

k moment does not necessarily exist for all basic distributions. 



63 



When it does, 



m k (Y » - (-D k ^E[e- sX ] 



s = -kA 



(1.14) 



For example, if X ~ Expon (1) , 



m k (Y) = 



k' 



(1 - kA) 



k+1 



if kA < 1; 



(1.15) 



otherwise the moment doesn't exist because the distribution's 
right tail is too long. 



(B) Symmetric Stretch-Producing Shapers 

Example: Z is a r.v. symmetrically distributed around 

2, 



zero; e.g., 



N(0,a ) . Here are some useful shapers: 
2 



(i) s(Z) = 1 + hZ~ , h > ; 

hz 2 
(ii) s(Z) = e (due to J. W. Tukey) . 



Again (i) and (ii) imply that Y = Zs(Z) resembles Z for 
small Z, but lengthens the tails of the distribution for 
large Z . 

Example 1 . Stretched log-Normal variables may be suggested 
for modeling repair or service times if data analysis indicates 
that the logarithms of observwd times are symmetric but not 
nearly Normal, having symmetrically too-long tails. Then it may 
be convenient to use the representation 



64 



Y = e 

(1.16) 
X = \x + cZs(Z) 

where u is the center (corresponds to mean) of the logged 
observations, and the sale constant c is the standard devi- 
ation (spread parameter) of the variable Zs(Z) , replacing o 
in the ordinary log-Normal formula. 
Moments . 

The moments of the shaped distribution are, for (i), 
representable in terms of those for a basic Z . 
(i) Y = Z(l + hZ 2 ) 



m 



1 (Y) = m (Z) + hm 3 (Z) , 

= (Z symmetrical around zero) . 



2 
m (Y) = m 2 (Z) + 2hm 4 (Z) + h m g (Z) . 



m 4 (Y) = m 4 (Z) + 4hm g (Z) + 6h 2 mg(Z) + 4h 3 m 10 (Z) 



4 
+ h m, 2 (Z) . 

? 2 

, ... „ hZ ; consider only Z ~ W(0,a ). Calculate 
li Y = Ze 



(1.17) 



as a preliminary 



/> z2 ,T" 2/ ° 2 _1_ - u _ Zha 2 )- 1 / 2 , h < 2a 2 (1.18) 

/2ttq 



65 



(1.19) 



Then differentiate repeatedly with respect to h to 
obtain the even moments (the odd moments equal zero) : 



2 

Var[Y] = m 2 (Y) = - 2 3/2 , h < l/4a" ; 

(1 - 4ha ) ' 



m 4 (Y > = ~ H 2,5/2 < h < i/ 80 ' • 

(1 - 8ha ) ' 

Hence kurtosis is 



m 4 (Y) 3(1 - 4ha 2 ) 3 
Y ? = — * = ^ o L 9 - 3 (1.20) 

(m 2 (Y)) Z (1 - 8ha^) *' z 



Tails are so extended that the kurtosis becomes very large in 
the case of (i) , and actually infinite for rather small values 
of h in (ii) ; the variance remains finite for slightly 
larger values. Nevertheless, the central part of the Y- 
distribution remains remarkably close to the Normal from when 
Z is itself Normal. 

Both forms (i) and (ii) can be induced to fit the 
inverse distribution (percent points) of the Student's t 
distribution fairly satisfactorily. 

(C) Left Tail Enhancement Shapers 

Example: X is a positive r.v., e.g. Exponential. 

(i) s(X) = 2£_ , i > o, a > 0; 

1 + aX 



66 



(ii) s(X) = 1 - e 



Ciii) s(X) = e~ a/X 



-aX' 



(1.21) 



This shaping function tends to concentrate probability near 
zero, leaving the distribution's shape unchanged for large X : 



Y = Xs(X) ~ < 



X 

(i) 
(ii) 



if X large 



V aX £+1 << X, X small 



(1.22) 



J 



Moments are generally impossible to find in useable form, 
although in the case of (ii) , £ = 1 , explicit results can be 
found in terms of Laplace transforms: suppose X ~ Expon(l) , 
then under (ii) 



E[Y] = E[X(1 - e aX ) ] = 1 - 



(1 + a) 



2 ' 



(1.23) 



the mean of Y is influenced very little when a is large, 
but small values of X are made even smaller; for the Expon(l) 
X the quantiles are related as follows: 



i \ r \ n aln(l-p) -, 

(p) = x(p) [1 - e ^ J / 



(1.24) 



so for a fixed p > large a forces y(p)/x(p) to one. 
However, for a fixed a and small p 



67 



*44~ ap 
x(p) - ^ 



(1.25) 



showing that low quantiles for X transform into even smaller 

values for Y . Furthermore, it is easy to show that if X,, » 

1 ( 1 ; n) 

is the minimum in a sample of n from X ~ Expon(l) , and 



-aX 



Y n \ = x ri x [1 - e 
(l;n) (l;n) L 



(l;n) 



then 



( 1 ; n) J 



a 



E[X 



(l;n) 



n + a 



as n -> oo 



(1.26) 



All of this reinforces the image of the present shaper as 
forcing small X-values to become smaller, leaving large values 
unchanged. Note that the result of the shaper (i) , £ = 1, 
can be explicitly inverted, giving the distribution 



F y (y) = F x 



"y(l + /1+4/ y) 



(1.27) 



from which the left-tail enhancement property shows itself 
explicitly: 



68 



f Y (y) 



= f 



x 



y(l + /1+4/gy 



h& 



/a /a /y 



f x (y) 



1 + y + C/a) 



y 2 - 4)y 



for y small 



for y large 



•(|) 



(1.28) 



Note that the hazard associated with the latter transformation 
is, when X is an Expon(l) r.v., 



h Y (y) = hi + y + v } 

/y + (4/a)y 



l//ay for y small 



for y large 



(1.29) 



Thus the density of an Exponential X remains, under this 
transformation, Exponential-like in the right tail but approaches 
infinity near the origin. A distribution with such behavior 
may well be useful for modeling failure data exhibiting early 
failures ("infant mortality.") 

(D) Right-Tail - Shortening Shaper 

Example: X is a positive r.v., e.g. an Exponential. To 
shorten or truncate the right tail, consider 



69 



(i) s(X) = 



1 + 3X 



9, 



3 > o, a > o 



(1.30) 



(ii) s(X) = e aX , 



a > 0/ £ > 



The shaper (i) with I = 1 is a right tail shortener or trunca- 
tor; it provides a monotonic transformation to Y <_ 3 ; small 
values of X are left nearly unchanged. The quantiles of Y 
are 

x(p) for p small 



^P> = 1 ^ X ( P ) - • 



(1.31) 



,-1 



for p large if x(p) >> 3 



-1 






Shaper (ii) is non-mono tonic: values of Y = Xs(X) increase 
to (ai) , decreasing thereafter. The inverse of (i) yields 



V^ = Vr^-s^ ' o i y i 3 



-i 



(1.32) 



with density 



f v ( ^ } = f x ( l - Bv } " 2 ; 

y XI 3y (1 _ &y) ^ 



(1.33) 



the latter approaches infinity at a rapid rate as y nears 
1/3 • 



70 



Shaper (i) with £ = 1/2 is a tail thinner; it can be 
inverted to give 



F v (y) = F Y [(By) 2 (^-^HIZE)2 ] . (1>34) 



with density 



f Y (Y) = f x C3y) 2 (il^tiVM) 2 ][i i (1+/1+(4/6 2 y) } (1 + i +(2/By) )3 

/l+(4/6 2 y 



f Y (y) for y small 



■x 



- { (1.35) 

2 2 2 

f (3 y )23 y for y large ; 



if X is Exponential, then the hazard rate of Y is initially 
flat, but eventually increases linearly with age, thus produc- 
ing a plausible wearout model for equipment or biological 
organisms. 

The above examples illustrate a few of the large number 
of possible ways in which the sculpturing idea can be used to 
extend the descriptive power of standard distributional models. 

Problems of fitting such models to actual data are 
currently being addressed, as are applications to simulation 
and time series modeling. Use of shaped exponentials to evaluate 
scheduling procedures has been initiated in collaboration with 
P. Collard at ULB. 



71 



(2) Response Times Under Processor Sharing 

Consider a simple model of a time-sharing computer sys- 
tem in which N terminals access a single computer (server) . 
Let repairman-model conditions prevail, but processor sharing 
governs the service order: if j jobs (0 <_ j <_ N) are at the 
execution stage each receives one-j — of a time unit of ser- 
vice per time unit. In other words, if it is the service rate 
(u is the expected service time under Markov-exponential 
assumptions) , then the unconditional probability that any des- 
ignated single job finishes in (t, t + s) is y(— ) + 0(A) . 

Under such conditions it is possible to derive backward- 
type equations to describe the response of waiting time, R , 
of a newly-initiated "tagged" job arriving from a previously 
idle terminal. In particular, consider 

m.(t) = E[R|W(R) = t , X(0) = j] , (2.1) 

where W(R) is the amount of work or processing time required 
of the server, and X(0) represents the number of jobs cur- 
rently at the processor when the tagged job first arrives 
(including that job). Thus m.(t) is the expected response 
time, conditional on need and accompaniment. Additionally, 
introduce r(j) as the fraction of a time quantum, D , actu- 
ally available for job processing when j jobs are present at 
the computer; r(j) represents one component of overhead, and 
may decrease as j increases. For short let A. represent 
the rate of new arrivals when j are present at the server; 



72 



under stated model conditions A. = A(N-j) . Let y. be the 
rate at which jobs accompanying the tagged job depart; under 

r j + 0(A) . These assump- 



stated conditions y. = y(j-l) 

tions, or those for a general birth and death process, lead 

by a backward-conditioning to this equation: 



m 



. Ct) = a + m .[ t - (^)a)[i - <A. + m .)a] + A.Am j + 1 



t ~ 



r(j) A 

j 



+ y . Am . , 



t - r <J> A 
j 



+ o(A) 



(2.2) 



.tract m. (t - ^4-^-) f: 
3 3 

let A -» to obtain 



Subt m. (t - " N y " ) from each side, divide by A , and 



^ dm.(t) 

— i = 1 - (A. + y ) m ( + ) + A. m..,(t) + y. m. . ( + ) (2.3) 

Ul - J J J J J -1 "- 1 - J j~ j - 



r(j) 

I 3 J 



This is a standard system of linear differential equations with 
constant coefficients that may be solved by standard methods. 

If r(j) is constant, and the repairman model assump- 
tions are fulfilled, then it has been shown by G. Latouche that 



e[r|w(r) = t] = c t , 



(2.4) 



i.e. is linear in t , with C depending upon A and y . See 
the article by Mitra [1981] for more detail concerning this problem, 



73 



( 3) Repairman Model in Random Environment 

Suppose we have m machines (this may mean electric 
power generators, or even remote computer terminals) that, when 
in use, fail independently (computer terminals: apply for 
processing time, or data) at rate A(t) , and, if failed are 
repaired at rate u(t) ; all processes are Markovian, given 
A(t) and u ( t) , t >_ . Now let A(t) and y(t) themselves 
be realizations of a finite-state Markov process that develops 
independently of the number N(t) of machines down for repair. 
If N(t) is the number down at t , and J(t) is the under- 
lying environmental state, then {N(t), J(t)} is a bivariate 
Markov process, and N(t) change is governed by the current 
level of the environment J(t) . The latter environment may 
refer to physical conditions such as heat, seismic shock, or to 
variations in repair effectiveness. In the case of computer 
terminals or communication nodes the environmental variations 
may be the result of changes in message transmissions or data 
demands under occasional crisis conditions. 

The paper by Gaver, Jacobs, and Latouche [19 81] presents 
a systematic mathematical analysis of the general birth-and- 
death process in random environments, including the above re- 
pairman model as a special case. Numerical illustrations are 
provided. Here we present a truncated version of the solution 
to the first-passage time problems, utilizing a recursive or 
"clawing-up" mode of thinking analogous to developments in 
Section ( ) of this account. Restrict discussion to just two 



7^ 



environmental levels, denoted by j = 1, 2 . Let the transi- 
tion rate from environmental state j = 2 -*■ 1 be a , and from 
j = 1 -+• 2 be 3 • Let ^ n (j) be the transition rate from 
(n,j) ■> (n+1, j) , and P n (j) be that from (n,j) ■*■ (n-1, j). 

Let, as before, U be the first-passage time from 
n to n + 1 , and put 



G n (dx;i,j) = P{U n e(dx), J (U n ) = j|N(0) = n, J(0) = i} , (3.1) 



with i, j = 1, 2 ; the L.-S. transform of this measure is 

called 

-sU 
G n (i,j;s) = E[e n , J(U n ) = j |N (0) = n, J(0) = i]; (3.2) 

this is the transform of the time to pass for the first time 
from a state in which n are down, the environment being in 
state i at some initial instant Ct = 0) , to n + 1 down, 
the environment being in state j . Simple considerations of 
cases that may arise during the first transition give, first 
for i = 1 , 

A (1) y (1) 2 

G n (1 '^ s) = d i riT £ j (1) + dnnr X G n -i (1 ' k ' s) G n (k ^ ;s) 

n n k — x 

(3.3) 



+ dTD G n < 2 '3;s) 



where j = 1, 2, the indicator function 



75 



if i = 

otherwise 



V j) = ■{ (3 - 4 > 



and the denominator 



d n (i) = A n (l) + y n Cl) + 3 + s . (3.5) 



Likewise, for i = 2 , 

X (2) u (2) 2 

G n (2 '^ s) = cTT2TV 2) + ^T2T I G n-l (2 ' k;s) G n (k '^ s) 
n J n k=l 

(3.6) 

+ dT2) G n (1 ^ ;s) 
n 

(the equivalent of these equations in Chu and Gaver [ ] 
accidentally incorrectly omits the final term on the rhs) . The 
above equations can in principle be solved recursively, 
beginning with 

A (1) a 



G o (1 '^ s) dTUT V 13 + dTOT G o (2 ' j;s) 

(3.7) 



'0 V " J ~0 



where 



G (2 '^ s) " dfuT £ j (2) + d^2T G (1 '^ s) 3 = 1.2. 



d Q (l) = A (l) + 3 + s . d Q (2) = X Q (2) + a + s . 



The first-passage time T R from N(0) = and J(0) = i (i=l,2) 
to n + 1 has transform 



76 



-sT 
P (i,j;s) = E[e n , J(T ) = j|N(0) = , J(0) = i] ; 



n - . n - 

in matrix notation, 

£ n - ie % • • • £ n 

and the L.-S. transform of the first-passage time to n + 1 

T 
in P £ , where Jl = (.1,1) , a column vector. 

Differentiation of expressions (3.2) and (3.5) produces 

recursive expressions for means, variances and higher moments. 

Programs have been written to evaluate these expressions 

numerically. 



77 



REFERENCES 



[1] AITCHISON J., and BROWN J. A. C. (1957) . The Log-Normal 

Distribution . Cambridge University Press, Cambridge, 
England. 

[2] BENDER E.A. (1978). An Introduction to Mathematical 
Modeling . John Wiley and Sons, Inc. , New York. 

[3] BURMAN D.Y. (19 79). An Analytic Approach to Diffusion 
Approximations in Queueing . Ph.D. Dissertation, 
NYU Courant Inst. 

[4] COX D.R., and SMITH, W.L. (1961). Queues . Methuen 
Monograph. John Wiley and Sons, Inc., New York. 

[5] COX D.R. (1962). Renewal Theory . Methuen Monograph. 
John Wiley and Sons, Inc., New York. 

[6] COX D.R. (1969) . Analysis of Binary Data . Chapman and 
Hall. London. 

[7] CHU B.B., and GAVER D.P. (.1977). Stochastic models for 

repairable redundant systems susceptible to common mode 
failure. S.I. A.M., Proc. of Int. Conf. on Nuclear 
Systems, Reliab. Eng. and Risk Assessment , Gatlinburg, 
Tenn. , 342-367. 

[8] CHUNG K.L. (1967). Markov Chains with Stationary Transition 
Probabilities. 2nd Ed. Springer-Verlag, New York. 

[9] CRAMER H. (1946). Mathematical Methods of Statistics . 
Princeton University Press, Princeton, NJ. 

[10] DANIELS H.E. (1954) . "Saddlepoint approximations in 
statistics". Ann. of Math. Statistics, Vol. 25, 
pp. 631-650. 

[11] FELLER W. (1951, 1966) . An Introduction to Probability 

Theory and Its Applications Vols. I and II . John Wiley 
and Sons, Inc., New York. 



78 



[12] GAVER D.P., and JACOB P. A. (1981). "On combinations of 
random loads". S.I. A.M. Jnl. Appl. Math. , Vol. 40, 
pp. 454-466. 

[13] GAVER D.P. and LEHOCZKY J. P. (1981). "Diffusion approxi- 
mations for the cooperative service of voice and data 
messages". J. Appl. Prob . , Vol. 18, pp. 660-671. 

[14] GAVER D.P., JACOBS P. A., and LATOUCHE G. (1981). "Finite 

birth-and-death models in randomly changing environments" 
Report Interne No. 121, Lab. D'Informat. Theor. , 
Univ. Libre de Brux. 

[15] GAVER D.P., and THOMPSON G.L. (1973). Programming and 
Probability Models in Operations Research . Brooks- 
Cole Publ. Co., Monterey, CA. 

[16] GAVER D.P., and SHEDLER G. (1973). "Approximate models 

for processor utilization in multiprogrammed computer 
systems". S.I. A.M. J. Comput. , Vol. 2, No. 3, 
pp. 183-192. 

[17] GAVER D.P. (1969) . "Highway delays resulting from flow- 
stopping incidents". J. App. Prob. , Vo. 6, pp. 137-153. 

[18] GAVER D. P. (1962) . "A waiting line with interrupted 
service, including priorities". J. of the Royal 
Stat. Soc , B, Vol. 25, pp. 73-90. 

[19] HARRISON J.M., and REIMAN M.I. (1981). "Reflected 

Brownian motion on an orthant" . Annals. Prob. , Vol. 9, 
pp. 302-308. 

[20] JACOBS P. A. (1980). "Heavy traffic results for single- 
server queues with dependent (EARMA) service and arrival 
times. Adv. Appl. Prob. , Vol. 12, pp. 517-529. 

[21] JACOBS P. A. (1978) . "A cyclic queueing network with 

dependent exponential service times". J. of Appl. Prob. , 
Vol. 15, pp. 573-589. 

[22] KARLIN S., and TAYLOR H.M. (.1975). A First Course in 
Stochastic Processes. Academic Press, London. 



79 



[23] KLEINROCK L. (1976) . Queueing Systems, Vols. I, II . 
John Wiley and Sons, Inc., New York. 

[24] KLINE, M.B., and ALMOG, R. [1930]. "Suitability o^ the 

loanormal distribution for repair times. Second Ins t. 
Conf. on Reliability and Maintainability Perros-Guirec- 
Tregastel, France, Sept. 1980. 

[25] KOBAYASHI H. (1981). This volume. 

[26] LEVY D. , and KAHN E. (1981). "Accuracy of the Edgeworth 
expansion of Loss-of-Load Probability calculations 
in small power systems". Summer Mtg., IEEE Power Eng. 
Soc. (IEEE Trans, on Power Apparatus and Systems, to 
appear) . 

[27] LEWIS P.A.W., and SHEDLER G.S. (1979). "Simulation of 

non-homogeneous Poisson processes by thinning". Naval 
Res. Log. Quart. , Vol. 26, No. 3, pp. 403-413. 

[28] LUGANNANI R. , and RICE S. (.1980). "Saddle point approxi- 
mation for the distribution of the sum of independent 
random variables". Adv. in Applied Prob. , Vol. 12, 
No. 2, pp. 475-490. 

[29] McNEIL D.R., and SCHACH S. (1973). "Central limit analogues 
for Markov population processes" (with discussion) . 
J. Roy. Stat. Soc. , Vol. 35, pp. 1-23. 

[30] McNEIL D.R. (1970). "Integral functionals of birth and 
death processes and related limiting distributions". 
Annals, of Statistics , Vo. 41, No. 2, pp. 480-485. 

[31] MITRA D. (1981). "Waiting time distributions from closed 
queueing network models of shared-processor systems". 
Proc. Eighth Internat. Symp. Comp. Performance Modelling 
Measurem. and Eval. Amsterdam, The Netherlands. 

[32] NEUTS M.F. (1981). Matrix-Geometric Solutions in Stochastic 
M odels - An Algorithmic Approach . The Johns Hopkins 
University Press, Baltimore Maryland. 

[33] NEUTS M.F., and MEIER K.S. (1981). On the use of phase type 
distributions in reliability modelling of systems with 
a small number of components. OR Spektrum, 2 , 227-234. 

80 



[34] NEWELL G.F. (1979). Approximate Behavior of Tandens Queues 
Springer-Verlag. Berlin, Heidelber, New York. 

[35] NEWELL G.F. (1971). Applications of Queueing Theory . 
Chapman and Hall. London, England. 

[36] PREGIBON D. (1981). "Logistic regression diagnostics". 
Ann, of Statistics , Vol. 9, pp. 705-724. 

[37] REIMAN M.I. (1982). Open queueing networks in heavy 
traffic" . To appear in Math, of Opns. Res. 

[38] TAKACS L. (1962). Introduction to the Theory of Queues 
Oxford University Press, New York, NY. 



DISTRIBUTION LIST 

NO. OF COPIES 

Library, Code 0142 4 

Naval Postgraduate School 
Monterey, CA 9 3940 

Dean of Research 1 

Code 012A 

Naval Postgraduate School 

Monterey, CA 9 39 40 

Library, Code 55 2 

Naval Postgraduate School 
Monterey, CA 93940 

Professor D. P. Gaver 14 8 

Code 55Gv 

Naval Postgraduate School 

Monterey, CA 93940 

Chief of NavalResearch 2 

Arlington, Virginia 22217 

Defense Technical Information Center 2 

ATTN: DTIC-DDR 
Cameron Station 
Alexandria, Virginia 22314 



82 



DUDLEY KNOX LIBRARY - RESEARCH REPORTS 



5 6853 01067370 



1)2022;