(* This notebook contains procedures that simulate various processes in the Gott-Caves dispute. The basic organization of this notebook is as follows.  In various comment cells, such as this one, we describe the simulations that the reader can perform.  Below each such cell, there is one or more procedure that is automatically initialized the first time that any command in this notebook is entered.  Examples of how these procedures are used follow the body of the procedure (in another cell). *)

(* This simulation consists of a sequence of busses, whose inter-arrival times are distributed according to one of the canned densities in Mathematica.  By this we mean that we imagine that a sequence of values from this distribution is created, and these values define the times between arrivals of busses at a certain bus stop.  This distribution is given the name cdist below, and can be chosen from the list given in th e statistical functions section of the Mathematica manual.  We imagine that the process has been running for a long time, so the process is near `equilibrium.' A would-be passenger appears at the bus stop, at some specific time.  The quantities that we will simulate are the time since the last bus arrived at the bus stop, the time until the next bus arrives at the bus stop, and the sum of these two times. We will call these quantities the age, the residual lifetime, and the lifetime, respectively, of the passenger's bus.  We denote the density of the inter-arrival times of the busses by f(t) and the cumulative distribution function of the inter-arrival time by F(t).  In renewal theory, it is shown that the density of both the age and the residual lifetime of the passenger's bus is equal to (1 - F(t))/m, where m is the mean of the inter-arrival distribution.  Also, it is shown that the density of the passenger's bus is given by tf(t)/m.  In the Gott-Caves dispute, it is important to understand the conditional distribution of the residual lifetime of the passenger's bus, given that the age of the bus equals some constant a.  By this we mean that we want to simulate the probability that the residual lifetime is at most y, given that the age equals a. In the Chance News article, we explain why this is equal to

F(y + a)/(1 - F(a)).  

Finally, Gott's main thesis is that if A represents the age of the passenger's bus, R represents the residual age of the passenger's bus, and c is a positive constant, then the probability that R exceeds cA equals 1/(1+c). *)

<< Statistics`ContinuousDistributions`

(* The above command is automatically executed when the user answers yes to the question asking whether the user wants all of the initialization cells entered.  This question is asked the first time any command is entered, once this notebook is open.  For example, to get this question to appear on the screen, one can put the cursor on the above command, and hit the `enter' key. *)

cdist  =  ExponentialDistribution[1]

PROGRAM:  Areabargraph
CALLING SEQUENCE:  Areabargraph[data, xmin, xmax, k]
PARAMETERS:
     data - a list of real numbers
     xmin, xmax - real numbers
     k - an integer
SYNOPSIS:
     -  This program divides the interval [xmin, xmax] on the x-axis into k subintervals
        of equal length.  Then, on each subinterval, a box is erected whose area is
        equal to the percentage of values in the list data which are in the subinterval.
        Important note:  only data values which fall in the user-defined interval will be
        included in	the bar graph.  If not all data values are included, a warning is
        displayed.  I f you wish to be sure of including all data values in the graph,
        choose xmin = Min[data] and xmax = Max[data].   (See also "Bargraph.")
RETURNED VALUES:
     -  none
LOCATION:
     File:  "Important Programs"

Clear[Areabargraph];
\(Areabargraph[data_, xmin_, xmax_, k_] := 
      Block[{sorteddata = {}, length = Length[data], 
          dx = \((xmax - xmin)\)/k, currentupperlim = xmin + dx, result = {}, 
          index = 1, counter, rectanglelist = {}, linelist = {}, 
          graphicslist = {}}, 
        sorteddata = 
          Sort[data]; 
If[\((\((sorteddata[\([1]\)] < 
                  xmin)\) || \((sorteddata[\([length]\)] > xmax)\))\), 
          Print["Note:  some data values lie outside the user-defined interval."]]; 
Block[{}, 
          For[i = 1, i ≤ k, \(i++\), 
            counter = 
              0; 
While[\((\((index ≤ 
                      length)\) && \((sorteddata[\([index]\)] ≤ 
                      currentupperlim)\))\), \(counter++\); 
\(index++\);]; 
result = 
              Append[result, counter]; 
currentupperlim += 
              dx;]; 
rectanglelist = 
            Table[Rectangle[{xmin + \((i - 1)\)*dx, 0}, {xmin + i*dx, 
                  result[\([i]\)]/\((length*dx)\)}], {i, 1, 
                k}]; 
rectanglelist = 
            Prepend[rectanglelist, 
              RGBColor[0, 0, 1]]; 
linelist = 
            Table[{Line[{{xmin + \((i - 1)\)*dx, 0}, {xmin + \((i - 1)\)*dx, 
                      result[\([i]\)]/\((length*dx)\)}}], 
                Line[{{xmin + i*dx, 0}, {xmin + i*dx, 
                      result[\([i]\)]/\((length*dx)\)}}], 
                Line[{{xmin + \((i - 1)\)*dx, 
                      result[\([i]\)]/\((length*dx)\)}, {xmin + i*dx, 
                      result[\([i]\)]/\((length*dx)\)}}]}, {i, 1, 
                k}]; 
linelist = 
            Flatten[linelist, 1]; 
graphicslist = 
            Union[rectanglelist, linelist]; 
Show[
            Graphics[{rectanglelist, linelist}, Axes -> True, 
              PlotRange -> All, Frame -> True]]]]\)

(* The following set of commands give basic examples of how one can use Mathematica to generate a set of random values from the distribution, put them in a histogram, and plot the theoretical density with the histogram. *)

Table[Random[cdist],  {100}]

(* It is sometimes helpful, for what follows, to know the maximum value in such a list.  The following command accomplishes this. *)

Max[%38]

(* The following command takes as input the list of values to be shown as a histogram (for example, the list above), the minimum and maximum values to be shown in the histogram, and the number of boxes into which the interval is to be cut. *)

Areabargraph[%38,  0,  8,  16]

(* The following command defines the theoretical density function.  This must be given by the user. *)

f[x_]  :=  E^(-x)

(* The following command plots the density function. *)

Plot[f[x],  {x,  0,  8}]

(* The following command shows the previous two plots on the same pair of axes. *)

Show[%45,  %47]

(* Now we wish to simulate the various quantities associated with the bus problem.  To make the process close to equilibrium, we allow the user to put in a time, arrivet, that the passenger arrives at the bus stop.  Presumably, the larger this time is (in relation to the mean of the interarrival times of the busses, the closer the process will be to equilibrium.  It should be pointed out that it is also true that the larger this time is, the longer this program takes to generate the required number of data values.  The number n, given by the user, is the desired size of the data set.  The output lists are called lastbuslist, nextbuslist, and passengerbustimeslist. *)

Bussimulation[cdist_,  n_,  
        arrivet_]  := 
Block[{lastbuslist  =  {}, 
              nextbuslist  =  {}, 
              passengerbustimeslist  =  {}, 
             timesubtotal, 
             currbustime
            }, 
            For[
          i  =  1,  
          i  ≤  
            n,  i++, 
                    timesubtotal  =  
            0; 
                     While[timesubtotal  <  
              arrivet, 
                                  currbustime  =  
              Random[cdist]; 
                                 timesubtotal  =  
              timesubtotal  +  
                currbustime;
                               ]; 
                    lastbuslist  =  
                       Append[
              lastbuslist,     arrivet  -  (timesubtotal  -  
                    currbustime)]; 
               nextbuslist  =  
                  Append[nextbuslist,  
              timesubtotal  -  
                arrivet]; 
                passengerbustimeslist  =  
                   Append[passengerbustimeslist,  
              currbustime];
                 ]; 
             Return[{lastbuslist,  
            nextbuslist,  
            passengerbustimeslist}]
            ]

(* The following procedure displays the output of the above procedure Bussimulation, along with the theoretical densities given by renewal theory.  The first two lists output by the above program should be distributed according to the inter-arrival density f(t).  The user should define this density BEFORE using this procedure.  The variable inputlists is just the set of three lists output by the above procedure Bussimulation.  The variables c1 and c2 are the endpoints of the interval on which the density is defined (c2 can equal Infinity) and c3 and c4 are the endpoints over which the lists will be graphed (both c3 and c4 must be finite).  The third list output by the above program is distributed according to the density t f(t)/m, where m is equal to the mean of the density (this is calculated by the procedure).

bussimulationplot[inputlists_, f_,  c1_,  c2_,  c3_,  
        c4_]  := 
Block[{lastbuslist  =  
            inputlists[[1]], 
             nextbuslist  =  
            inputlists[[2]], 
             passengerbuslist  =  
            inputlists[[3]], 
             m, 
           }, 
           m  =  
          Integrate[t f[t],  {t,  c1,  c2}]; 
Show[
          Areabargraph[lastbuslist, c3,  c4, 
            20],  
          Plot[(1  - 
                  Integrate[f[s],  {s, c1, t}])/m,  {t,  c3, 
              c4}]
        ]; 
Show[
          Areabargraph[nextbuslist, c3,  c4,  
            20],  
          Plot[(1  - 
                  Integrate[f[s],  {s, c1, t}])/m,  {t,  c3, 
              c4}]
        ]; 
Show[Areabargraph[passengerbuslist, c3,  c4,  20],  
          Plot[t*f[t]/m,  {t,  c3,  
              c4}]];
              ]

(* The first distribution we will try is the exponential density with parameter lambda = 1.  *)

cdist  =  ExponentialDistribution[1]

f[t_]  :=  E^(-t)

(* The following set of calculations are for an arrival time of 50. *)

Bussimulation[cdist,  10000,  50];

bussimulationplot[%40,  f,  0,  Infinity,  0,  10]

(*The next distribution we will try is the uniform density with on the interval [0, 1]. *)

f[t_]  :=  1

cdist  =  UniformDistribution[0,  1]

Bussimulation[cdist,  10000,  500];

bussimulationplot[%80,  f,  0,  1,  0, 1]

(* Now we wish to simulate the statement made in the Chance News review concerning the distribution of the remaining time (until the next bus), given that the time since the last bus is known.  It is difficult to simulate this exactly for continuous distributions, since we will never be able to generate data for which the time since the last bus equals the given value.  The way out is to sort into boxes the values of the time since the last bus, and use the midpoint value of each box to represent the time since the last bus.  The next procedure does exactly this. It assumes that some data has already been produced by the preceding program.  The user declares a value of a and the radius r of an interval, centered at a, which defines which data points are relevant.  The tradeoffs here are obvious; if r is small, then the data more accurately approximates the assumption that the time since the last bus equals a, but there are fewer data points in the sample that fit in the interval. The output is a list of the values of the remaining time (the time until the next bus).  The user must input the two lists (which are the first two output by the preceding program), a value of a, and a value of r.  This procedure will then print the number of data points satisfying the given constraint, and return the list described earlier. *)

\(nextbusgivenage[lastbuslist_, \ nextbuslist_, \ a_, \ 
        r_]\  := 
Block[{consolidatelist\  = \ {lastbuslist, \ 
              nextbuslist}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ subsetlist\  = \ {}
\ \ \ \ \ \ \ \ \ \ \ }, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ subsetlist\  = \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \(Transpose[
              Select[Transpose[
                  consolidatelist], \ \((\((#[\([1]\)]\  ≥ \ 
                          a - \ r)\) && \((#[\([1]\)]\  ≤ \ 
                          a + r)\))\)\  &]]\)[\([2]\)]; 
Print["The size of the output list is ", \ 
          Length[subsetlist]]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ Return[subsetlist]
\ \ \ \ \ \ \ \ \ \ \ \ \ \ ]\)

nextbusgivenage[%69[[1]],  %69[[2]],  1,  .1];

Max[%]

Areabargraph[%111,  0,  5,  20]

Show[%113,  %47]

(* The final quantity that we will simulate is the probability that R > cA, and plot the results as a function of c.  On the same pair of axes, the theoretical value 1/(1+c) for this probability is also plotted.  The values c1 and c2 are defined by the user, and represent the endpoints of the interval over which the theoretical probability will be plotted.  The values c1 = .1 and c2 = 10 usually work well. On a color monitor, the data curve is in red and the theoretical curve is in blue. *)

residualversusage[lastbuslist_,  nextbuslist_,  c1_, 
        c2_]  := 
Block[{n  =  
            Length[lastbuslist], 
            i,  
             roveralist  =  {}, 
             fraclist  =  {}, 
             }, 
              roveralist  =  
          Sort[nextbuslist/
              lastbuslist]; 
              For[i  =  1,  
          i  ≤  
            n,  i++, 
                      fraclist  =  
              Append[fraclist,  {roveralist[[i]], (n - i)/
                    n}];
                  ]; 
Show[
          ListPlot[fraclist,  Axes  ->  True, 
            AxesOrigin  ->  {0, 0},  PlotJoined  ->  True,  
            DisplayFunction  ->  Identity,  
            PlotStyle  ->  {RGBColor[1, 0, 
                  0]}], 
    Plot[
            1/(c + 1),  {c,  c1,  c2}, Axes  ->  True,  
            AxesOrigin  ->  {0, 0}, 
            DisplayFunction  ->   Identity,  
            PlotStyle  ->  {RGBColor[0, 0, 
                  1]}], 
DisplayFunction  ->  $DisplayFunction]; 
Return[
          fraclist]   
             ]

residualversusage[%69[[1]],  %69[[2]],   .1,  
        10];

(* The following procedure picks n real numbers uniformly from the interval [0, m].  These real numbers can be viewed as starting times for busses.  Next, the lifetimes of these busses are picked from a distribution with density given by cdist.  Finally, we collect all of the lengths of the busses that are alive at time m/2.  We display these with an areabargraph.  The user must define, before the procedure below is called, the theoretical density f(t) for the interarrival times of the busses.  The user should also give the interval [c1, c2] on which the density is defined (the value Infinity is allowed). *)

Clear[specialbusses];
specialbusses[n_,  m_,  cdist_,  f_,  c1_,  
        c2_]  := 
Block[{buslist  =  {}, 
              randlife, 
             randstart,  
             i, 
             mean 
           }, 
             For[
          i  =  1,  
          i  ≤  
            n,  i++, 
                     randlife  =  
            Random[cdist]; 
                     randstart  =  
            Random[Real,  {0,  
                m}]; 
                     If[((randstart  <  m/2) && (randstart  +  randlife  >  
                    m/2)), 
                           buslist  =  
              Append[buslist,  
                randlife]
                        ]
                   ]; 
mean  =  
          Integrate[t*f[t],  {t,  c1,  c2}]; 
Show[
          Areabargraph[buslist, c1,  Min[c2,  10],  20],  
          Plot[t*f[t]/mean,  {t,  c1,  
              Min[10, c2]}]];
  ]

cdist  =  ExponentialDistribution[1]

f[t_]  :=  E^(-t)

specialbusses[10000,  10,  cdist,  f,  0,  Infinity]

f[t_]  :=  1

cdist  =  UniformDistribution[0,  1]

specialbusses[10000,  10,  cdist,  f,  0,  1]