This example shows how to solve a statistical simulation and display problem. Imagine that telephone calls arrive at a bank of n ports and each call is assigned to the first open port found by scanning, in order, ports 1, 2, .... We assume n is large enough so that no call will be denied a port. Suppose that the length of the ith call, l_i, follows a lognormal distribution, so that, for all i, the length, l_i, of the ith call has the probability density function

f(x) = exp(-log(x/m)2/(2s2))/(xs(2p)1/2),

and

P(li £ h)= ó õ h 0  f(x) dx.

The average length of a call is then m·exp(s²/2), and we will choose m = 10 and s = .5. The call lengths l₁, l₂, ... are all independent. Also, suppose calls arrive at varying times with exponentially-distributed inter-arrival times. Let A_0j be the time between the arrivals of the jth and (j+1)st calls. Then we assume that the probability that there are s or fewer minutes between the arrival of the jth call and the arrival of the (j+1)st call is P(A_0j £ s) = 1-e^-as, where a is, say, 1 minute.

It is of interest to simulate this process and to graph the simulated use of the various ports as follows. In the graph below, the time of arrival of calls are marked on the t-axis, and each use of port i is shown by a solid line, starting at the indicated time and extending for the duration of the call.

A graph of this form is determined by a list of triples, one for each arriving call; so that for call i, we have (t_i, l_i, p_i) where t_i is the time at which call i arrives, l_i is the length of call i, and p_i is the port to which call i is assigned. Then we need only draw the line segments (t_i,p_i) to (t_i+l_i,p_i) for i = 1, 2, ¼, and, of course, mark the t-axis at the points (t_i, 0) to indicate the times at which the calls arrive.

Given data of this form, we can, for example, easily estimate the percentage of time port i is busy, or estimate the number of ports needed to service 95% of the incoming calls.

Let A_ij denote the interarrival time, A_ij, between the jth and (j+1)st calls assigned to port i; A_ij is a random variable. An analysis of interest is to state the distribution of A_ij in terms of i, j, and the distribution functions of the length and interarrival time of incoming calls. It often suffices to determine the limiting form for the distribution of A_ij as j®¥.

In the case where the times of arrival of incoming calls forms a Poisson process as assumed above, we can explicitly state the desired distributions.

Let D_ij be the non-busy waiting time between the end of the jth call on port i and the arrival of the next call assigned to port i. Let L ~ l₁, so L is lognormally-distributed, as is l₁. Note the inter-arrival time A_ij ~ L+D_ij for i ³ 1, and also D_1j ~ A_0j.

Let D_i=lim_j®¥D_ij. Indeed D_ij ~ D_i for j > 1. Then, in the limit, port i is busy the fraction f_i of the time, where f_i = E(L/(L+D_i)).

Now, when port i is free, an arriving call will be assigned to port i when ports 1, 2, ..., i-1 are all busy. So if a call arrives during an interval of time which satisfies this condition, it will be assigned to port i. Thus, in the limit, D_i is the waiting time until the occurrence of the event that ports 1, 2, ..., i-1 are busy and a call arrives, and this can be interpreted as the waiting time until the next event of a Poisson process with density f₁f₂¼f_i-1/a, where 1/a is the density of incoming calls. Thus D_i is exponentially-distributed with mean a/(f₁f₂¼f_i-1), for i > 1, and D₁ ~ A_0j.

Moreover, in the limit, the number of busy ports, N, at a fixed time is a Poisson random variable with mean u = E(L)/a; that is P(N = k) = u^k e^-u/k!.

Note this model is of interest where "calls'' are sessions with a time-sharing computer, people arriving for service at a system of queues which are scanned linearly, and other isomorphic situations. Let the vector (V₁,V₂,¼,V_n) denote the state of the n ports at the current time, t, so that V_i = 0 if port i is not busy, and V_i = the length of time remaining for the call using port i otherwise. Now, given that the next new call of length l arrives in s minutes from the current time t, we can update the state vector V to be valid at time t+s-e, where e is very small, by changing each V_i so that V_i ¬max(V_i-s,0). Then the state vector V is updated to be valid at time t+s+e by setting V_j to l, where j = X(V), and X(V) is the least index in {1,2,¼,n} such that V_X(V) = 0. Let us denote the updating process used to make V the current state vector for time t+s-e, so that we may write V ¬ u_s(V) to denote updating V to time t+s-e.

Now our simulation can be programmed as follows. Note a total of d calls are simulated.

begin
    t₀ ¬ 0; j ¬ 1; V¬ 0; a:t_j ¬ random-exponential+t_j-1;
    l_j ¬ random-lognormal;
   V ¬ u_tj-tj-1(V); p_j ¬ X(V); V_pj ¬ l_j; j ¬ j+1; if j £ d then goto a;
end.

Upon completion, we have the three vectors, t, l, and p needed to draw the required graph.

We may use the computer program MLAB to carry out the simulation sketched above and then draw the generated graph of port-usage over time. The MLAB system was originally developed at the National Institutes of Health and includes curve-fitting, differential equations, statistics and graphics as some of its major capabilities. MLAB is a tool for researchers in science and engineering. MLAB is a tool for solving continuous simulation and parameter-estimation problems such as chemical kinetics and compartmental models as well as discrete simulation problems such as we are considering here. A general do-file (script) to do the simulation and draw the graph of port-usage is given below, together with a sample graph generated by using it.

"calls.do" 
type enter n: the number of calls"; 
n = kread();

type enter a random number seed"; 
r = kread(); r = ran(r);

type enter np: number of telephone ports"; 
np = kread();

v = 0^^np; "(np ports. Port i is busy for v[i] minutes)";

"generate t vector in order: partial sums of samples of an exponential 
distribution with mean 1"; 
t = -(log on (ran on 0^^n)); 
t = 0 psum(t);

"generate l vector: lognormal with mean 10*exp(-(1/8)).";
l = normran on 0^^n; l = 10*(exp on .5*L);

"compute the p vector";
function pv(i) = if i > np then np else if v[i] = 0 then i else pv(i+1);
function u(w) = if w > s then w-s else 0; 
p[n] = 0; 
for j = 1:n do  \
   {i = minrow(v); if v[i] not= 0 then {np=np+1; i=np;}; 
   p[j] = i; v[i] = l[j]; 
   s = (if j = n then 0 else t[j+1])-t[j];
   v = u on v}

"given t, l, and p, draw the graph."; 
m1 = t &'p; m2 = (t+l) &'p;
draw m1 line none pointtype vbar color red in w1 
draw mesh(m1,m2) line alternate color green in w1 
draw t '0 line none pointtype vbar in w1 
draw 0 '(1:np) line none label 1:np offset (-.02,-.01) color yellow in w1 
no imagebox in w1 
view

Here is a graph of a simulation run with 100 calls produced by the above do-file.

This simulation can be run in a variety of forms. For example, we may assume that the mean of the exponential distribution for the inter-arrival times of incoming calls varies smoothly with the time of day, thus mimicking a heavier load during working hours, with peaks at, say, 10:00am and 3:00pm. This can be done by assuming the mean inter-arrival time, E(A_0j), is an appropriate doubly-peaked function of time.