Patent Publication Number: US-2007124123-A1

Title: Computer system for simulating a physical system

Description:
FIELD OF THE INVENTION  
      The present invention relates to a computer system for simulating a physical system.  
     BACKGROUND  
      In a clinical trial, an experimenter administers a fixed dose of a drug to several patients and monitors the concentration of the drug in the blood stream as a function of time. The clinical trial can be simulated as a model using a computer system to determine whether the concentration of the drug exceeds a therapeutic level, but falls below a harmful level.  
      The model may depend on several factors such as the rate at which the drug passes from the stomach to the blood stream and the rate at which it is excreted. These rates vary from person to person. Moreover, concentration depends on the volume of distribution, which in turn depends upon body weight. The model uses probability distribution functions of rates and weights. However, due to the possible high number of distributions and their complexity, plentiful computational resources are needed to model the trial. Otherwise it can take a long time to simulate the trial. In fact, this problem is common in simulating physical systems.  
      The present invention seeks to provide an improved computer system in which a physical system can be simulated more quickly and/or with less processing power, while maintaining accuracy.  
     SUMMARY OF THE INVENTION  
      According to a first aspect of the present invention there is provided a computer system for simulating a physical system comprising a first module configured to convert a probability distribution function of a physical property into a discrete function having at least one segment and a second module configured to perform an exact transformation of each segment.  
      Each segment may have a non-zero width and an approximated probability density. The probability density may be constant across at least one of the at least one segment. The probability density may vary across at least one of the segments. The first module may be configured to ensure that the sum of the areas of the at least one segment equals one. The first module may be configured to convert the probability distribution function into a piecewise-defined function having only one segment. The first module may be configured to convert the probability distribution function into a piecewise-defined function having more than one segment and to set a boundary between adjacent segments dependent upon rate of change of the probability distribution function.  
      The second module may be configured to add results of each transformation to provide a new probability distribution function. The second module may be configured to perform a unary transformation of each segment. The second module may be configured to perform a binary transformation of each segment with each segment of another discrete probability distribution function, such as convolution of each segment with each segment of another discrete probability distribution function and/or addition, subtraction, multiplication or division integral transformation of each segment with each segment of the other discrete function.  
      The binary transformation may produce, for each pair of segments, a function having first and second regions, wherein, in the first region, the probability density increases from zero to a non-zero value and, in the second region, the probability falls from a non-zero value to zero. The function may have a third region, between the first and second regions. The functions arising from each pair of segments may be added to produce a new probability distribution function. The binary transformation may produce, for each pair of segments, a triangular function having first and second regions, wherein, in the first region, the probability increases linearly from zero to a non-zero value and, in the second region, the function falls linearly from the non-zero value to zero.  
      The first module may be configured to convert the new probability distribution function into a new discrete function having at least one segment. The second module may be configured to perform another exact transformation of each segment of the new discrete function. The other exact transformation may differ from the exact transformation.  
      The probability distribution function may be of an amount of a drug or a time interval in a processing system.  
      According to a second aspect of the present invention there is provided a method of simulating a physical system, the method comprising converting a probability distribution function of a physical property into a discrete function having at least one segment and performing an exact transformation of each segment.  
      According to a third aspect of the present invention there is provided a computer program product comprising a computer-readable medium storing a computer program for causing a computer to convert a probability distribution function of a physical property into a discrete function having at least one segment and to perform an exact transformation of each segment. 
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS  
      Embodiments of the present invention will now be described, by way of example, with reference to the accompanying drawings in which:  
       FIG. 1  are illustrative plots of concentration of a drug in the blood stream as a function of time;  
       FIG. 2  is a schematic block diagram of a computer system for simulating a physical system in accordance with certain embodiments of the present invention;  
       FIG. 3  is a process flow diagram of a simulating process performed by the computer system shown in  FIG. 2 ;  
       FIGS. 4   a  to  4   d  illustrate, for a simple example, results of performing steps in the simulating process shown in  FIG. 3 ;  
       FIGS. 5   a  and  5   b  illustrate, for another example, results of the simulating process shown in  FIG. 3 ;  
       FIG. 6  illustrates a unary transformation;  
       FIG. 7  illustrates a binary transformation;  
       FIG. 8  is a plot of number of requests in a processing system which can be modelled in accordance with certain embodiments of the present invention; and  
       FIG. 9  is a plot of waiting time of the requests shown in  FIG. 8 . 
    
    
     DETAILED DESCRIPTION OF EMBODIMENTS  
      Drug companies usually formulate a drug such that when it is administered to a patient its concentration, for example in the blood stream, exceeds a level at which it provides a therapeutic effect, but falls short of a harmful level. However, the effect of the drug will vary from person to person, for example due to differences in their body weights. Thus, drug companies usually formulate the drug to have a dose such that the concentration falls within these levels when taken by most of the population. The size of dose is usually found as part of a clinical trial.  
      For example, referring to  FIG. 1 , if a fixed dose of a drug is administered to three different patients, then the concentration of the drug in their respective blood streams could vary with time as shown in plots  1   1 ,  1   2 ,  1   3 . For example, this may be due to differences in body weights and/or other factors.  
      Before conducting any trial, an experimenter may wish to ascertain ranges of doses which are safe to use. Furthermore, during the trial, the experimenter may wish to check that the results of the trial are as expected. For this reason, the experimenter simulates the concentration of the drug in the blood stream based on an understanding of how the body processes the drug.  
      Referring to  FIG. 2 , a computer system  2  for simulating a physical system, such as concentration of a drug in the blood stream, in accordance with certain embodiments of the present invention is shown. The computer system  2  may be distributed and, thus, may include a client  3  and server  4  connected via a network  5 . The client  3  may be a relatively low-specification personal computer, whereas the server  4  can be a high performance computer. Alternatively, the computer system  2  may be implemented in a single computer.  
      The computer system  2  includes a module  6  for supplying at least one input probability distribution function  7 , a module  8  for supplying at least one rule  9  for processing the at least one input probability distribution function  7  and a module  10  for outputting and post-processing at least one output probability distribution function  11 . These modules  6 ,  8 ,  10  may be implemented on the client  3 .  
      Each probability distribution function  7  defines a variable or a parameter of a physical system, such as body weight.  
      The probability distribution function(s)  7  and rules  9  may be prepared manually by the user and/or selected from a pre-stored database.  
      The computer system  2  further includes a first module  12  for converting the at least one input probability distribution function  7 , if not already discrete (herein also referred to as “piecewise-defined”), into corresponding discrete or piecewise-defined probability distribution function(s)  13 , a module  14  for performing an exact probabilistic transformation of the discrete probability distribution function(s)  14  so as to provide at least one transformed probability distribution function  15  (herein also referred to as the “result”) and a second module for converting the result  15  into at least one discrete result  17 .  
      The first discretising module  12  can take a probability distribution function  7 , such as f(x), and divide at least part of the function  7  over an interval (x 1 →x 2 ) into several segments over a plurality of intervals (x 1 →x 1,1 , x 1,1 →x 1,2 . . . , and x 1,n →x 2 ). The discretising module  13  approximates the probability distribution function  7  in each or, if there is only one, the segment, and may ensure that the sum of the areas of the segments equals one. However, the discretising module  12  can simply approximate the probability distribution function  7  over the interval, such that there is a single segment.  
      The discretising module  12  divides the probability distribution function  7  according to the rate of change of the probability distribution function  7  so as to provide higher resolution when the probability distribution function  7  changes rapidly and lower resolution when the probability distribution function  7  changes slowly. However, the discretising module  12  can divide the probability distribution function  7  into equidistant intervals or equiprobable intervals. The number of segments can be chosen to provide a given degree of accuracy. Typically, the number of segments is of the order 1 to 10,000. However, the number of segments can be higher.  
      The discretising module  12  approximates the probability distribution function  7  in each interval as a uniform probability distribution, i.e. with constant probability. Thus, the discrete, approximate probability distribution function  13  resembles a histogram or bar chart. However, other probability distributions, such as non-uniform linear or quadratic probability distribution functions, may be used provided they can be exactly (e.g. analytically) transformed.  
      The transformation module  14  may perform unary transformations, such as square root or cosine, of a function or binary transformations, such addition, subtraction, multiplication or division integral transforms, of the probability distribution function(s)  7 . Other transformations are possible. Furthermore, one or more unary and one or more binary transformations may be combined, for example sequentially.  
      The transformation module  14  is able to perform (or perform more easily) an exact, for example analytical, probabilistic transformation because the discretising module  12  has converted each probability distribution function  7  into discrete, uniformly valued segments.  
      The second discretising module  16  (hereinafter referred to as the “repartitioning module  16 ”) takes what may be a continuous result  15  and converts it into a discrete result  17  which can be fed back into transformation module  14  and/or output.  
      The computer system  2  may be provided with an optional database  18  for storing probability distribution functions  7 ,  13 , rules  9  and results  15 ,  17 .  
      The computer system  2  includes at least one processor (not shown), memory (not shown) storing software for providing the modules  6 ,  8 ,  10 ,  12 ,  14 ,  16 , device drivers (not shown), input devices (not shown), such as a keyboard and/or mouse, output devices (not shown), such as a display, and storage devices (not shown), such as hard disk.  
      Referring to  FIG. 3 , a process flow diagram of a process of a method of simulating a physical system in accordance with certain embodiments of the present invention is shown.  
      Referring also to  FIG. 4   a , the process will be described using simple probability distribution functions  7   1 ,  7   2 .  
      The probability distribution functions  7  and rules  9  are specified by a user and supplied by the probability distribution function-supplying module  6  ( FIG. 2 ) and rule-supplying module  8  ( FIG. 2 ) to the discretising module  12  ( FIG. 2 ) and transformation modules  14  ( FIG. 2 ) respectively (steps S 301  &amp; S 302 ).  
      A first probability distribution function  71  is a uniform function over a first interval  20 . A second probability distribution function  72  is an exponential function.  
      The discretising module  12  ( FIG. 2 ) divides an interval  21  of the second probability distribution function  7   2  into first and second intervals  21   1 ,  21   2  and approximates the probability distribution function  7   2  in each interval  21   1 ,  21   2  to provide a discrete function having first and second segments  22   1 ,  22   2  (step S 303 ).  
      In this example, the rule  9  is to convolve the first and second probability distribution functions  7   1 ,  7   2 , that is: 
 
 h ( x )=ƒ( x )⊕ g ( x )  (1) 
 
      Referring to  FIG. 4   b , the transformation module  14  ( FIG. 2 ) applies the rule  9  to produce a result  15  ( FIG. 2 ), which in this example involves transformation module  14  ( FIG. 2 ) convolving probability distribution function  7   1  with the segments  22   1 ,  22   2  producing partial solutions  23   1 ,  23   2  which, when added, produce a cumulative result  23  (step S 304 ).  
      Referring to  FIG. 4   c , a comparison of the cumulative result  23  and a full, conventional analytical result  24  is shown. It will be appreciated that this comparison is for the purposes of illustration only and that the computer system  2  does not usually carry out a full, conventional analytical solution not least because the solution is impossible or onerous to compute.  
      The repartition module  16  ( FIG. 2 ) may convert the result  15  ( FIG. 2 ), in this case cumulative result  23 , into a discretised result  17  ( FIG. 2 ) (step S 305 ).  
      Referring to  FIG. 4   d , the repartition module  16  ( FIG. 2 ) may divide the cumulative result  23  into intervals  25   1 ,  25   2 ,  25   3  and approximate the cumulative result  23  in these intervals  25   1 ,  25   2 ,  25   3  as approximations  26   1 ,  26   2 ,  26   3 .  
      The process may continue to process the result  15  (step S 306  &amp; S 303 ). For example, the transformation module  14  ( FIG. 2 ) may perform further transformations on approximations  26   1 ,  26   2 ,  26   3 .  
      The results may be output to the output module  10  ( FIG. 2 ) for displaying and post-processing by the user (step S 306 ). Post-processing may include determining quantiles, determining mean and variance and/or determining area under the curve.  
      As mentioned earlier, the computer system  2  ( FIG. 2 ) can simulate the concentration of the drug in the blood stream.  
      The concentration, C, is the amount of drug in the blood stream, B, divided by the so-called “volume of distribution”, V, i.e.  
             C   =     B   V             (   2   )             
 
      For a single individual, the volume is given: 
 
V i =k i W i   (3) 
 
 where W i  is the body weight of the individual and k i  is a person-specific constant. 
 
      Thus, the concentration for the individual i is:  
               C   i     =     B       k   i     ·     W   i                 (   4   )             
 
      For an individual, B, k i , W i  are single numbers and so it is easy to decide whether C i  falls between a therapeutic level T and a harmful level H.  
      For a population, however, it is harder to decide whether the concentration C i  falls with these limits. Within a population, body weight has a probability density PDF(W) and constant also has a probability density PDF(k). Thus, for a given dose, not all individuals will have the same concentration C i . Even if the mean concentration &lt;C i &gt; lies within limits, the concentration for some individuals may fall below the therapeutic level T or exceed the harmful level H. By calculating PDF(C), the experimenter can assess the efficacy of the clinical trial.  
      To simulate the concentration of the drug in the blood stream, the probability distribution function module  6  ( FIG. 2 ) supplies the following probability distribution functions  7  ( FIG. 2 ) to the discretising module  12  ( FIG. 2 ): 
 
V e =normal (0, 0.18)  (5a) 
 
w=normal (70, 20)  (5b) 
 
d=normal (100, 0.1)  (5c) 
 
b=uniform (0, 0.1)  (5d) 
 
e=uniform (0, 0.1)  (5e) 
 
k 1 =normal (1, 0.1)  (5 ) 
 
k 2 =normal (1, 0.1)  (5 g) 
 
 where V e  is the random variation in the volume of distribution, w is body weight, d is the amount of drug in the stomach, b is the amount of drug in the blood stream, e is the amount of drug in the kidneys (i.e. excreted), k 1  is the decay constant for the drug passing from the stomach to the blood stream and k 2  is the decay constant for the drug passing from the blood stream to the kidneys. 
 
      In equations 5a to 5f above, two distributions are used, namely normal(μ,σ 2 ) which is the normal distribution P(x) in a variate x having mean μ and standard deviation σ 2  and defined by:  
               P   ⁡     (   x   )       =       1     σ   ⁢       2   ⁢   π           ⁢     ⅇ         -       (     x   -   μ     )     2       /   2     ⁢     σ   2                   (   6   )             
 
 and uniform(a, b) which is the uniform distribution function defined by:  
               P   ⁡     (   x   )       =     {         0           for   ⁢           ⁢   x     &lt;   a               1   /     (     b   -   a     )               for   ⁢           ⁢   a     ≤   x   ≤   b             0           for   ⁢           ⁢   x     &gt;   b                     (   7   )             
 
      The discretising module  12  ( FIG. 2 ) converts the normal probability distribution functions  7  into discrete functions having, for example, three segments (not shown).  
      The rule-supplying module  8  ( FIG. 2 ) supplies a function, namely: 
 
V=0.21We V     e     (8) 
 
 where 0.21 is phenomenological constant. 
 
      The rule-supplying module  8  ( FIG. 2 ) also supplies differential equations, namely:  
                 ⅆ   d       ⅆ   t       =       k   1     ⁢   d             (     9   ⁢   a     )                   ⅆ   b       ⅆ   t       =         k   1     ⁢   d     -       k   2     ⁢   b               (     9   ⁢   b     )                   ⅆ   e       ⅆ   t       =       k   2     ⁢   b             (     9   ⁢   c     )             
 
      The transformation module  14  ( FIG. 2 ) applies the rules found in equations 8 and 9a to 9c, to the discretised probability distribution functions based on equations 5a to 5e.  
      Different output may be selected. For example, PDF(d), PDF(b) or PDF(e) can be selected at different times t, for example t=0.2, 0.4, 0.6, 0.8, 1.0 and 1.2.  
      Referring to  FIG. 5   a , a screen shot  27  of an output provided through the output module  10  ( FIG. 2 ) is shown. The screen shot  27  shows a plot  28  of PDF(b) at t=1.0.  
      Referring to  FIG. 5   b , a screen shot  28  of another output provided through the output module  10  ( FIG. 2 ) is shown. The screen shot  29  shows plots  30   1 ,  30   2 ,  30   3 ,  30   4 ,  30   5  of C=b/V with time t. The plots  30   1 ,  30   2 ,  30   3 ,  30   4 ,  30   5  include not only a plot  30   3  of medium concentration, but also quantile plots.  
      As explained earlier, the transformation module  14  ( FIG. 2 ) is able to perform an exact probabilistic transformation because the discretising module  12  ( FIG. 2 ) has converted each probability distribution function  7  into discrete, uniformly valued segments.  
      Some examples of exact unary and binary transformations of discrete probability distribution functions will now be described in more detail.  
      Referring to  FIG. 6 , an example of a discrete probability distribution function  31  is shown. The discrete probability distribution function  31  may correspond to a discrete probability distribution function  13  ( FIG. 2 ) produced by the discretisation module  12  ( FIG. 2 ).  
      The probability distribution function  31  has three segments  31   1 ,  31   2 ,  31   3  each having respective variate starting and finishing values of (x 1 , x 2 ), (x 2 , x 3 ) and (x 3 , x 4 ) and probability distribution values a 1 , a 2  and a 3 .  
      A unary transformation of the discrete probability distribution function  31  is carried out by transforming each segment separately  31   1 ,  31   2 ,  31   3 . A square root transformation will be used to illustrate a unary transformation:  
      The transformation module  14  ( FIG. 2 ) takes the square root of the starting variate value  32 , (x 1 ) of the first segment  31 , to produce a new starting variate value x 1   ′ =√{square root over (x 1 )} and the square root of the finishing variate value  32   2  (x 2 ) of the first segment  31   1  to produce a new finishing variate value x 2   ′=√{square root over (x 2 )}. Between the new starting and finishing variate values x   1   ′ , x 2   ′  the new segment  31   1   ′  is given a renormalized probability distribution value a 1 . The same process is applied to the second segment  31   2  and so on.  
      Referring to  FIG. 7 , further examples of discrete probability distribution functions  33 ,  34  are shown. The discrete probability distribution functions  33 ,  34  may correspond to a discrete probability distribution functions  13  ( FIG. 2 ) produced by the discretisation module  12  ( FIG. 2 ).  
      A first probability distribution function  33  has one segment  33   1  having variate starting and finishing values of (X A , X B ), in this case (1, 2), and probability distribution value c, in this case 1. A second probability distribution function  34  has one segment  34   1  having variate starting and finishing values of (X C , X D ), in this case (3, 5), and probability distribution value j, in this case 0.5.  
      A binary transformation of the discrete probability distribution functions  33 ,  34  is carried out to produce a transformation probability distribution function  35 . An addition integral transformation of the discrete probability distribution functions  33 ,  34  will be used to illustrate a binary transformation. The transformation module  14  ( FIG. 2 ) determines a starting variate value  36   1  , of the transformed probability distribution function  35  by adding the start variate values of the functions  33 ,  34 , i.e. X A +X C . The transformed probability distribution function  35  rises up a first slope  37   1  to a first shoulder variate value  36   2 . The transformation module  14  ( FIG. 2 ) may determine the first shoulder variate value  36   2  using, for example: 
 
X A +X C +min(X B −X A , X D −X C ) 
 
      The transformed probability distribution function  35  may have a plateau  38 , although it need not necessarily do so. The transformation module  14  ( FIG. 2 ) may determine a second shoulder variate value  36   3  using, for example: 
 
X B +X D −min(X B −X A ,X D −X C ) 
 
      The transformed probability distribution function  35  drops down a second slope  37   2  to a finishing variate value  36   4 . The transformation module  14  ( FIG. 2 ) determines the finishing variate value  36   4  by adding the finishing variate values of the functions  33 ,  34 , i.e. X B +X D .  
      If the first and second discrete probability distribution functions  33 ,  34  have more than two segments, then this process is repeated for each pair of segments. For example, this involves convolving a first segment from the first probability distribution function  33  with a first segment from the second probability distribution function  34 , the first segment from the first probability distribution function  33  with a second segment from the second probability distribution function  34  and so on until the first segment from the first probability distribution function  33  with a last segment from the second probability distribution function  34  then a second segment from the first probability distribution function  33 , with the first segment from the second probability distribution function  34  and so on.  
      For an addition integral transformation, the slopes  37   1 ,  37   2  are linear. For other forms of integral transformation, the slopes  37   1 ,  37   2  mall have other shapes, such as logarithmic. The plateau  38  need not have a constant value, but may have, for example in the case of multiplication integral transformation, a logarithmic shape.  
      The system and process hereinbefore described may be used to simulate any system which is driven by randomness, such as queuing or manufacturing processes.  
      For example, the system may be a computer system (hereinafter referred to as a “processing system”), in which requests are processed using a First-In-First-Out (FIFO) arrangement.  
      A request j arrives at a system time t j . If the system is busy, then the request waits for a time w j . The system services the request in time b j , after which a result leaves the system at time r j . The time span between arrivals j−1 and j may be denoted a j  and the time span between departures j−1 and j may be denoted d j . The time span for which the system is idle before the arrival of request j may be denoted i j .  
      Individual events and distributions of event are listed in Table 1 below:  
                   TABLE 1                       Individual events   Distribution of events                                                t i     arrival time   —   —       w i     waiting time   W   distribution of waiting times       b i     service time   B   distribution of service times       r i     departure time   —   —       a i     inter-arrival time   A   distribution of inter-arrival times       d i     inter-departure time   D   distribution of inter-departure times       i i     idle time   I   distribution of idle times                  
 
      Referring to  FIG. 8 , a plot  39  of the number of requests waiting to be processed by a processing system is shown. The plot shows periods  40   1 ,  40   2  when the system is busy and a period  41  when the system is idle.  
      Referring to  FIG. 9 , a plot  42  of waiting time is shown. The plot  42  includes lines  42   1 ,  42   2 ,  42   3 ,  42   4  representing a time a request will wait in the system before departure. A new arrival is represented by a jump of magnitude b j .  
      If the system has been idle (as for j=1 and j=4), then the remaining service time for the request immediately starts to decrease. However, if the system has been busy (as for j=2 and j=3), then a newly-arrived request is subject to a waiting time which equals the sum of the service times of previous requests.  
      Thus, the following relations hold: 
 
 r   j−1   =t   j−1   +w   j−1   +b   j−1   (10) 
 
 r   j   =t   j   +w   j   +b   j   (11) 
 
 and therefore  
               d   j     =         r   j     -     r     j   -   1         =         (       t   j     -     t     j   -   1         )       ︸     a   j         +     (       w   j     -     w     j   -   1         )     +     (       b   j     -     b     j   -   1         )                 (   12   )                 w   j     =     max   ⁡     (     0   ,       w     j   -   1       +     b     j   -   1       -     a   j         )               (   13   )             
 
 and a so-called “breathing space” can be defined as: 
 
 c   j−1   =b   j−1   −a   j   (14) 
 
 such that equation  13  can be re-written as 
 
 w   j =max(0 ,w   j−1   +c   j−1 )  (13) 
 
 and 
 
 d   j   =b   j −min(0 , w   j−1   +c   j−1 )  (15) 
 
 i   j =−min(0 ,w   j−1   +b   j−1   −a   j )  (16) 
 
      For a given sequence of inter-arrival times and service times, a sequence of waiting, inter-departure times and idle times can be determined.  
      Equations 13′, 15 and 16 can be translated into operations for probability distribution functions, namely: 
 
 C ( t )= B ( t )⊕ A (− t )  (17) 
 
 W ( t )=π( W ( t )⊕ C ( t ))  (18) 
 
 D ( t )= B ( t )⊕ I ( t )  (19) 
 
 I ( t )=+E,vos π( W ( t )⊕ C ( t ))  (20) 
 
 where π and  π  are sweeping operations for maximum and minimum respectively and where equation 19 can be presented as:  
               D   ⁡     (   t   )       =       ∫     -   ∞       +   ∞       ⁢       B   ⁡     (   t   )       ⁢     I   ⁡     (     t   -   τ     )       ⁢     ⅆ   τ                 (     19   ′     )             
 
      Equations 18 and 20 can be presented in a similar format.  
      The π operation involves integrating the probability distribution function W(t)⊕C(t) for negative variate values (for example, over the range −∞→0) so as to find a weighting α and replacing this part of the probability distribution function with a delta function at zero variate value, the delta function having pre-weighting α.  
      The  π  operation involves integrating the probability distribution function W(t)⊕C(t) for positive variate values so as to find a weighting β and replacing this part of the probability distribution function with a delta function at zero variate value, the delta function having pre-weighting β. The negative part of the probability distribution function is flipped by reflecting the function about the y-axis.  
      The waiting time probability distribution function Wis determined self consistently. This may involve guessing an initial waiting time probability distribution function W, such as there is no waiting, and iteratively calculating equation 18 until the result converges. For example, this may involve determining whether two results are the same, i.e. W i =W i−1 , a or any difference is small, for example based on chi-squared test. Once the waiting time probability distribution function W has been determined, the inter-departure times and idle times probability distribution functions D, I can be determined.  
      Thus, for a given inter-arrival time probability distribution function A and a service time probability distribution function B, probability distribution functions W, D, I for waiting, inter-departure times and idle times can be determined.  
      The computer system  2  ( FIG. 2 ) can thus be used to simulate the processing system and obtain probability distributions functions for waiting, inter-departure times and idle times.  
      The system and process hereinbefore described may be used to simulate any system which is driven by randomness. The process has the advantage that they can be used to simulate a wider class of stochastic models than existing analytical methods and are more accurate and faster at computing models than using Monte Carlo methods.