\(\(Get["VisualStats`"]; 
Needs["Statistics`Master`"]; 

    Needs["Graphics`Master`"]; 
$DefaultFont\  = \ {"Courier", \ 8}; 
    
\)\)

MATHEMATICA IN ACTION for Issue 8-1 of Mathematica in Education and Research

Statistics in the Classroom and the Courtroom

The field of statistics has long been known for innovative use of graphics. Mathematica’s breadth allows us to create graphics suited to specific needs, from the education of students to presentations to a jury.

by Karl Heiner and Stan Wagon

.  Introduction

The abstract nature of some of the fundamental concepts in statistics causes difficulty for many students. Experimentation and case studies are pedagogical tools that facilitate the study of these concepts. Both provide concrete examples to which to tie the concepts and both require that students remain “on task” in order to complete the assignment. It is well known by experts in education that on-task behavior is a required component of effective learning. This is especially critical for students who have grown up in a world of sound bites and commercial TV. In this column we describe some approaches to teaching statistical concepts that show how one might use the power and versatility of Mathematica to develop functions that allow students to experiment and analyze case studies. Moreover, the examples are based on actual data sets and cases, thus providing a solid link from abstract theorems and distributions to the real world.

.  Hypothesis Testing

.1. An Employment Discrimination Case

In the early 1990s New York State was experiencing a serious budget crisis. As a result, Governor Cuomo instructed each division of the state government to lay off approximately 10% of its employees. So the Management Administration Unit of the Office of Mental Health dismissed 11 of its 127 employees. The average age of the terminated employees — call it avg — was greater than than that of those not laid off. Those terminated filed an age discrimination suit that was argued in the Federal District Court in Albany. The ages of the individuals in each group are given below.

\(Needs["Statistics`Master`"]; 
Needs["Graphics`Master`"]; 

    Needs["VisualStats`"]; \)

\(kept = {29, 34, 36, 44, 36, 47, 33, 33, 45, 45, 36, 
		26, 29, 37, 
          29, 30, 38, 38, 57, 40, 47, 41, 42, 45, 48, 
		33, 62, 43, 44, 
          32, 40, 32, 25, 47, 36, 31, 24, 50, 26, 
		42, 45, 32, 46, 28, 
          42, 53, 45, 57, 26, 62, 46, 49, 47, 
		31, 40, 42, 28, 45, 30, 
          53, 27, 48, 29, 40, 48, 47, 42, 
		44, 36, 66, 40, 40, 47, 38, 
          29, 31, 37, 39, 47, 42, 58, 
		60, 52, 48, 43, 47, 37, 39, 21, 
          41, 31, 39, 39, 53, 41, 
		25, 48, 51, 41, 23, 48, 48, 62, 66, 
          44, 48, 35, 45, 37, 
		46, 34, 40, 33, 32, 40, 30} // N; 

    dismissed = {42, 39, 37, 57, 33, 42, 60, 39, 56, 44, 28} // N; 

    all = Join[kept, dismissed]; \)

In such cases there is an initial burden on the plaintiff to show that there was disparate impact on the protected class; this is usually accomplished by a statistical test of the hypothesis of no disparate impact. The following demonstration was used by K.H. to introduce the concept of hypothesis testing to the jury and to demonstrate that obtaining an average age at least as large as avg when a random drawing of eleven from the combined list was not particularly unusual. First, we note the size of the groups and their average age.

Map[Length, {kept, dismissed}]

{avgkept, avg} = Map[Mean, {kept, dismissed}]

Using the RandomKSubset function from the VisualStats package we can easily get the ages of eleven randomly chosen individuals.

RandomKSubset[all, 11]

We can generate as many of these random layoff lists as we wish, say five, and find the mean age for each one.

Mean /@ Table[RandomKSubset[all,  11], {5}]

By generating 100 random layoff lists we get an idea of the likelihood of randomly picking a list with an average age equal to or greater than that of the plaintiffs (43.3636). This time 21 of the 100 means exceeded 43.3636, which is not rare enough for the jury to conclude disparate impact. Inasmuch as the plaintiffs’ lawyers failed to meet this burden of proof, much of their remaining evidence was irrelevant.

SeedRandom[1] ;
means = Map[Mean, Table[RandomKSubset[all,  11], {100}]];

bc = BinCounts[means, {32, 50}]

Count[means, x_ /; x > avg]

\(\(BarChart[bc, PlotRange -> {{0, 18}, {0, \ Max[bc] 1.1}}, 
	
      BarLabels -> Range[33, 50], BarStyle -> {Gray}, 
      Frame -> True, 
	FrameTicks -> None, Axes -> True, 
	
      Epilog -> {Text["avg", {avg - 31.8, 17}, {\(-1\), 0}], 
          AbsoluteThickness[4], \ Line[{{avg - 32, 0}, {avg - 32, 17}}]}, 
	
      Ticks -> {\(({#, # + 32}\ &)\)\ /@\ Complement[Range[2, 17], {3}], 
          Range[2, 18, 2]}, AspectRatio -> 1/3, 
      AxesOrigin -> {0, 0}]; \)\)

Actually, the argument that the average age of those on the layoff list was greater than the average age of those not on the list was itself irrelevant, because by law the protected class is those over the age of 40. The usual statistical test employed by the plaintiff’s attorney to argue that the protected class has been disparately impacted is the Fisher exact test, which is based on the hypergeometric distribution; see section 3.1.2.

.. The Fisher Exact Test and the Hypergeometric Distribution

In cases such as the preceding, ordinarily the plaintiff will construct a 2×2 table displaying the numbers of individuals included on and excluded from the layoff list by class status (protected vs. not protected) and attempt to argue that, because of their disproportionate representation on the layoff list., the protected class was treated unfairly. Then, using the hypergeometric distribution, one calculates the probability of randomly observing the table in question or one more extreme (i.e., more discriminatory), given that the marginal counts are held constant (e.g., 55 less than 40 years of age and 72 older, 11 listed for layoff, while 116 were not). If the sum of these probabilities is small, we reject the hypothesis of no disparate impact. This is called Fisher’s exact test. By observing the table for our age discrimination case (see below), it is clear that a lesser proportion of the protected class was listed for layoff than that proportion of the unprotected class. Here there is no need to calculate probabilities, since it is clear that there is a greater proportion of people over 40 in the group not laid off than there is in the laid off group. More about the hypergeometric distribution later.

Code for Table (do not print, adjust font size as appropriate)

\(fontsize = 9; 

    DisplayForm[
      StyleBox[GridBox[
          Prepend[{
		
              Prepend[Append[
                  b1 = BinCounts[dismissed, {39.5 - 100, 39.5 + 100, 100}], 
                  Plus@@b1], "Dismissed"], 
		
              Prepend[Append[
                  b2 = BinCounts[kept, {39.5 - 100, 39.5 + 100, 100}], 
                  Plus@@b2], "Kept"], 
		
              Prepend[{
                  b1[[1]] + 
                    b2[[1]], 
                  b1[[2]] + 
                    b2[[2]], ""}, 
                "Total"]}, 
            
{"", "<40", "≥40", "Total"}], 
\ 
          GridFrame -> 1, RowLines -> {1, 0, 1}, 
          ColumnLines -> 1], 
FontFamily -> Times, 
        FontSize -> fontsize, Background -> GrayLevel[0.9]]]\)

.. A Randomization Test Lab Activity

Scheaffer et al [Scheaffer 1996] describe a number of discrimination cases that lend themselves to analysis using randomization tests. For example, in a 1977 race discrimination case, Dendy vs. Washington Hospital Center, four of nine African-American nurses passed an examination while all 26 white nurses passed. In such cases it may be argued that employment tests are discriminatory if they adversely impact a protected class. To be successful, the plaintiff endeavors to shift the burden of proof to the defendant by using statistical evidence to reject the null hypothesis that the impact on the protected class is the same as that on the unprotected class. If the plaintiff is successful in establishing this disparate impact, then the defendant must prove that the use of the criterion is a business necessity.

An ordinary deck of playing cards can help us develop a reference distribution for assessing the likelihood of the observed event: four of nine African-American nurses passed an examination while all 26 white nurses passed. If we discard 17 black cards from the deck and shuffle it thoroughly, we can deal five cards representing one realization of this examination situation where success on the exam is unrelated to race. The dealt cards correspond to individuals who fail. If we repeat this experiment many times we can assess the likelihood that all five failures are from the nine nurses belonging to the protected class.

In teaching the concepts of hypothesis testing and reference distributions (or working in court) it is useful to have the students (or jurors) actually perform the experiment with a deck of cards; but such physical experiments waste time and are distracting.  We can use the VisualStats package function Deal. Evaluate the following ten times (or even twenty) and you will almost surely see no cases where five black cards turn up. The RemovedCards option sets the removed cards in the suits clubs, diamonds, hearts, spades (in that order, using {1, 2, …, 13} for the cards.) The cards themselves were designed using Bézier curves as discussed in [Wagon 1995].

Deal[5, ColorCount -> True, 
      RemovedCards -> {Range[13], {}, {}, Range[4]}];

\(\(Show[
      GraphicsArray[
        Map[First, 
		
          Table[Deal[5, DisplayFunction -> Identity, 
				
              RemovedCards -> {Range[13], {}, {}, Range[4]}], {5}], 
          2]]]; \)\)

Code for Figure

\(SeedRandom[1]; 
lightgray\  = \ GrayLevel[0.4]; 


    Show[GraphicsArray[
        Map[First, \ 
            Table[
			
              Deal[5, 
			\ 
                DisplayFunction\  -> 
							Identity, 
	
                RemovedCards -> {Range[13], {}, {}, Range[4]}, \ 
                FontSize -> 10], {5}], 2] /. \ 
          RGBColor[1, 0, 0] -> lightgray]]; \)

Of course, we could have avoided experimentation by applying our knowledge of the hypergeometric distribution from our study of discrete probability distributions. The probability that four of nine African-American nurses passed an examination while all 26 white nurses passed is given by the hypergeometric distribution (see section 3.1.2). Such a calculation shows that the chance of the observed event is about 1 in 2600, which provides good evidence that the test is discriminatory. More about the hypergeometric distribution in a moment.

.  Some Probability Models

The understanding of the axioms of probability and some special probability distributions is essential to the understanding of hypothesis testing, and of course these probability models have important applications in their own right. In this section, we demonstrate the use of some VisualStats experimentation functions that we ask students to use to understand and reinforce how these models might naturally arise.

..  Models Arising from Sampling

Here we describe a few experiments to help the student understand the similarities and differences between the binomial and hypergeometric random variables. These random variables arise naturally from sampling.  Although these two distributions are a result of attribute sampling, the notions about sampling carry over to variable sampling.

...  Attribute Sampling With and Without Replacement

Some people distinguish between attribute sampling, where one samples individuals or sampling units from some population and then observes the presence or absence of some attribute (e.g., over 40 or not, black or not, success or failure, defective or not), and variable sampling, where one takes measurements of a more continuous nature (e.g., height, weight, dollars, calories). When the population from which the sample is drawn is finite, we distiguish between sampling with replacement and sampling without replacement. To illustrate, consider an urn containing r red balls and b black balls. We will choose one ball at random, record its color, and then sample again after having replaced the original ball (or not). A more general replacement strategy replaces the selected ball together with c balls of the color selected and d balls of the other color. For the hypergeometric distribution, which assumes sampling without replacement, c=-1. The binomial distribution requires sampling with replacement, or c=0. In both cases, d=0.

The VisualStats function Urn provides a convenient means for performing this experiment. For more on urn models see [Feller, 1968]. Urn[{r, b}, {c, d}, n] returns images of n generations of an urn that starts with r red balls and b black balls. At each generation, a ball is chosen and c balls of the same color of the chosen one and d of the other color are added to the urn (c can equal -1).

For example, let’s return to the case of 26 white nurses and 9 African-American nurses who attempted an exam that produced five failures. If we assume that success on the exam is inpedendent of race and that 30 jobs are available (i.e., five must fail), then the process can be described using an urn model and sampling without replacement. Of particular importance is the fact that the probability of observing a red ball changes from trial to trial.

The code that follows generates a 5-generation movie of an urn evolution (see Figure ).

Urn[{26, 9}, {-1 , 0}, 5];

Code for Figure

\(\(SeedRandom[1]; 

    images\  = \ 
      Block[{$DisplayFunction\  = \ Identity}, 
		
        \((\ uu = 
            Urn[{26, 9}, {\(-1\), 0}, 6, \ BallColors -> {White, Black}, 

              LabelBackground -> White, \ 
              DefaultFont\  -> \ {"Times", \ 8}, StepSize -> 1])\); 
		
        tt = Table[
            uu[\([i]\)]\  /. \ \ 
              \((PlotLabel -> z_)\)\  :> \ 
				
                \((PlotLabel\  -> 
                    DisplayForm@
						\ \ 
                      \((StyleBox[
                          GridBox[
                            Map[MakeBoxes[#]\ &, 
										
                              Partition[
                                Cases[z, \ 
                                  StringForm[x_, y__]\  :> \ 
                                    StringForm[
                                      StringReplace[x, \ 
                                        "Red" -> "White"], y], \ 
                                  ∞], 2], {2}], \ GridFrame -> 1], 
                          \ FontSize -> \((9)\), \ FontFamily\  -> \ Times])
                        \))\), \ {i, \ 6}]\ ]; 
\)\)

\(Off[Show::shx]; 

    Show[GraphicsArray[
        \ {
			{images[\([1]\)], \ images[\([2]\)]}, 

          	{images[\([3]\)], \ images[\([4]\)]}, 

          		{images[\([5]\)], \ images[\([6]\)]}}]]; \)

...  Hypergeometric Distribution

Suppose we had a population of size N consisting of k successes, and consequently N-k failures.  If we draw a sample of size n without replacement, the probability of exactly x successes and n-x failures is

(((({{k}, {x}}))(({{N-k}, {n-x}})))/((({{N}, {n}})))).

If a sample of size five is drawn without replacement from a population of size 35 consisting of 9 black balls and 26 red ones, the probability of five blacks, and consequently 0 reds, is

N[((Binomial[9, 5] Binomial[26, 0])/Binomial)[35, 5]]

Or, we may use functions from Mathematica’s standard package Statistics`DiscreteDistributions`. The hypergeometric distribution is the probability function returning the probability of exactly x successes in a sample of size n drawn without replacement from a population of size N comprised of k successes and N-k failures.  The function PDF[HypergeometricDistribution[n, k, N], x] from the DiscreteDistributions package provides this probability mass.

N[PDF[HypergeometricDistribution[5, 9, 35], 5]]

As expected from our investigation into the nurse-exam simulation, this is an extremely rare event.

...  The Binomial Experiment

We can return to the urn model and now require that the selected ball be replaced. Unlike sampling without replacement, in this situation the probability of obtaining a red ball remains the same from trial to trial. By specifying a fixed number of draws or trials, n, and counting the number of successes (red balls) observed, x, we have the binomial experiment. The probability of a success on each trial, p, is merely the relative frequency of red balls in the urn; p can of course be any number from zero to one, 0≤p≤1. Sampling with replacement in the urn model can be seen to be equivalent to observing the number of heads in a sequence of coin flips.

The function CoinFlips[n, 1] used next performs a useful experiment. Students can experiment with it by changing the number of tosses or the relative frequency (use the HeadProbability option) of red balls.

Realistic examples are useful. Suppose that one of the products produced by a computer manufacturer is very new and requires a component, a particular computer chip, that is very difficult to produce in large quantities. However, a supplier of the chip has been found, but the supplier will only assure that 85% of the chips function reliably. In chip manufacturing this is known as the process yield. Defective chips will be detected at one stage of the manufacturing process and the cost of replacing the defective chips in terms of time and materials has been factored into published prices. It is critical that shipments of the chips contain no more than 15% defective, so procurement has recommended a lot acceptance sampling procedure that will test some number — say n — of chips from each lot. In order to eliminate the cost of returned shipments, the testing will take place at the supplier’s facility.

If the manufacturer directs the supplier to randomly choose seven chips for testing and refuses to accept a lot if the seven tested chips include more than one defective chip, what proportion of the time will procurement refuse a shipment with 85% good chips? Here we may ask whether the sampling should be with or without replacement, which distribution is appropriate, and when does the choice of distribution matter.

Execute the next cell for an example of what might happen when the results of the tests are recorded. On the inspector’s recordkeeping form there is box representing each chip to be inspected.  If a chip is defective, he darkens the corresponding box on his form. Try this a few times. We may change the number of chips inspected and/or the proportion defective.

SeedRandom[1]; 
 CoinFlips[7, 1, MovieOnly -> True, HeadProbability -> 0.15]
        ;

This problem is analogous to flipping seven coins. Experiments similar to a coin-tossing experiment occur frequently.  And the CoinFlips function can be used in several ways to study such binomial experiments. By setting the MovieOnly option to False, the default, and by using the second argument to indicate the number of times to repeat the binomial experiment, we introduce the notion of a sampling distribution and the central limit theorem. Notice that to the right of each realization of the binomial experiment is printed the numbers of heads, which are in turn tabulated in the histogram on the right. The HeadProbability option can be used when the probability of a head is different than one half.

Through experimentation, the student should discover that as the number of coins gets large, the histogram displaying the number of heads resembles the Gaussian distribution; this is the DeMoivre–LaPlace theorem, a special case of the central limit theorem. The student will also notice that rate of convergence is faster when the probability of a head is closer to one half.

CoinFlips[7, 15];

Code for Figure

SeedRandom[1]; 
    CoinFlips[7, 15, TextStyle -> FontColor -> Black];

.. Distribution Associated with Waiting Times

Experimenting with randomly occurring events in time allows us to explore the Poisson and exponential distributions. The Poisson random variable is a discrete random variable; the exponential random variable is continuous. We discuss them together because in some applications they are intimately related. In such cases, each of these random variables is defined on the same sample space, the sample space resulting from a Poisson experiment. These two random variables are not always naturally connected and may be applied in a variety of situations, but here they will be introduced together as random variables evolving from the Poisson process.

... The Poisson Process

The Poisson process describes randomly occurring events in a continuum, such as time or space (of 1, 2, or 3 dimensions). For the moment we will use time as an example of a continuum in which random events may occur. In the demonstration that follows, events are randomly generated in time. You can think of these events as customers arriving at a bank, or cars arriving at a toll booth. Note that we may view this experiment as random events in time and are interested in how many events occur during some time period, or the amount of ellapsed time between arrivals.

The input cell below generates a series of graphics cells that form an animation demonstrating the Poisson process, a process where events appear randomly in time. The images may be interpreted as people walking up to a teller’s window. The number of arrivals in any time period is a Poisson random variable. The distribution of the time intervals between arrivals will follow the exponential distribution.

WaitingTimes[20, Movie -> True]

... The Poisson Random Variable and Its Probability Mass Function

In the Poisson process where events occur randomly in time, the number of events occuring in some fixed time interval follows the Poisson probability distribution. That is, if we fix the time interval over which we intend to observe the process, the random variable X that counts the number of events occuring in this interval is the discrete random variable called the Poisson random variable. The probability mass function for this random variable is f(x)=(e^(-λ)) (λ^x)/ x!, where λ is the average number of events occurring in a time interval of this length.

Suppose that customers arrive randomly at a bank at a rate of 30 per hour. Thus if we are interested in observing the process for one hour, λ is 30. On the other hand if we wish to observe for one half hour, λ = 15. What is the probability of exactly 10 arrivals in a one half hour?

N[((15^10)/((E^15) 10!))]

In its DiscreteDistributions package, Mathematica provides a function for the Poisson distribution. We can view its equation for the probability mass (or density) function (pdf), by evaluating the cell below.

PDF[PoissonDistribution[λ],x]

And we may use these functions to obtain the probability of exactly x events. In the preceding example we wanted the probability of exactly 10 arrivals.

N[PDF[PoissonDistribution[15], 10]]

We can view the shape of this probability mass function by constructing a bar chart to display the individual probabilities.

SeedRandom[1]

BarChart[Table[{PDF[PoissonDistribution[15], x], x}, {x, 0, 30}], 
      PlotRange -> All, BarStyle -> {GrayLevel[0.5]}];

This bar chart may be somewhat misleading in the sense that the probability of any noninteger is zero, but the bars may give the illusion that there is some nonzero probability for values of X close to, but not equal to the integers. We attempt to represent these point masses by lines with nearly no width.

SeedRandom[1]

\(λ = 15; n = 30; 

    Plot[0, {x, \(-1\), n + 1}, 
      Epilog -> 
        Table[Arrow[{x, 0}, {x, PDF[PoissonDistribution[λ], x]}, 
            HeadShape -> 
              Point[{x, PDF[PoissonDistribution[λ], x]}]], {x, 0, 
            n}], PlotRange 
        -> {{\(-1\), n + 1}, {0, 
            PDF[NormalDistribution[λ, \@λ], λ]}}, 
      AxesLabel -> {"x", "f(x)"}, 
      AxesOrigin -> {\(-1\), 0}]; \)

To obtain the probability of no more than 10 arrivals, we must sum the probabilities of 0, 1, …, 10.

Sum[N[PDF[PoissonDistribution[15], x]], {x, 0, 10}]

Of course, what we were asking for in the cell above was the cumulative distribution function or cdf of 10, F(10).

N[CDF[PoissonDistribution[15], 10]]

The mean and the variance of the Poisson distribution are both λ. If λ gets large, the shape of the Poisson distribution gets more symmetric and we get the impression that we could use the normal distribution to approximate it, using the appropriate correction for continuity.

Clear[λ];
{Mean[PoissonDistribution[λ]],
Variance[PoissonDistribution[λ]]}

The Poisson random variable is a discrete random variable since there is only a countable infinity of possible values of the count of the randomly occuring events in time. If, on the other hand, we are interested in the continuous random variable that measures the distance between randomly occurring events in the Poisson process, the exponential distribution may be derived.

... The Exponential Random Variable and Its Probability Density Function

When observing a Poisson process, we can, instead of counting events in some fixed time interval, record the time between successive events. The probability distribution of the continuous random variable T, the time between successive events, frequently follows the exponential probability distribution, whose probability density function, {t,f(t)} is defined by the relationship f(t)=(1/β) (e^(t/β)) for t>0. Here β represents the average time between arrivals of these randomly occurring events in time, and therefore is the reciprocal of λ, the mean of the Poisson distribution. This fact may become more obvious if we consider our customer arrivals at the bank that averaged 30 per hour. The expected time between arrivals is thus 2 minutes: β = 1/λ.

We can plot this probability density function and notice that there is an uncountably infinite number of possible values of our random variable T. In order to find the probability that T is between two values a and b, we find the area under the pdf between a and b.

λ = 2; Plot[
      PDF[ExponentialDistribution[λ], t], {t, 0, (4/λ)}];

The expected value and variance of the the exponential distribution are β and (β^2), respectively.

Clear[λ];
{Mean[ExponentialDistribution[λ]],
Variance[ExponentialDistribution[λ]]}

We may wish to convince ourselves that (e^(-t/β))/β is truely a pdf (nonnegative and the total area under this curve is 1).

Integrate[
      PDF[ExponentialDistribution[(1/β)], t], {t, 0, ∞}, 
      Assumptions -> Re[β] > 0]

We can find the probability that T<c, for any positive number c. We merely find the area under the pdf below the point c on the t-axis.

\(Clear[β]; 

    ∫\_0\%c PDF[ExponentialDistribution[1\/β], t] 
        \[DifferentialD]t\)

Recall that we call this the cumulative distribution function or cdf of c, or F(c).

CDF[ExponentialDistribution[1/β], c]

So, we discover that we do not need to use calculus to find areas under the exponential pdf; we merely make use of the fact that F(c)=1-(e^(-c/β)). The next code generates an exponential pdf and shades the area below the value 1.5. The shaded area is F(1.5)=1-(e^(-1.5/λ)). The area above any point c is just (e^(-c/β)).

\(λ = 2; 

    Show[{FilledPlot[PDF[ExponentialDistribution[λ], t], {t, 0, 1.5}, 
          PlotRange -> {{0, 4\/λ}, {0, 2}}, \ 
          DisplayFunction -> Identity], 

        Plot[PDF[ExponentialDistribution[λ], t], {t, 0, 4/λ}, 
          PlotRange -> {{0, 4/λ}, {0, 2}}, \ 
          DisplayFunction\  -> Identity]}, 

      DisplayFunction -> $DisplayFunction]; \)

... Queueing

One very useful application of these two probability distributions is in the area of the theory of waiting lines, or queues.

The code that follows produces a graphic representation of a simulation of a waiting line situation. In this simulation, 10 customers arrive at a server, say a drive-up teller. The time to serve these customers is 1.5 times the interarrival time. For each of the 10 bars produced in this graphic, there is a black segment at the left end of the bar which is the same for each customer and represents the time to get from the street to the end of the line. On the right end of each bar is a black (or green) segment depicting the service time for that customer. Between the drive-up segment and the service segment, there is the potential of a gray (or red) waiting-time segment. In this example notice that the mean service time is 1.5 times the mean interarrival time, so a waiting line tends to develop. The vertical line indicates when the window closes. Below the horizontal bars is a graph that depicts the length of the waiting line (the queue) at any time. Change the service time (the BetaValue option) and watch the impact on the queueing.

WaitingTimes[10, BetaValue -> 1.5, Servers -> 1];

Code for Figure

\(SeedRandom[1]; 

    ww = WaitingTimes[10, BetaValue -> 1.5, Servers -> 1, \ \ 
        DisplayFunction -> Identity]; 


    Show[ww\  /. 
        \ {RGBColor[1, 0, 0]\  -> \ GrayLevel[0.5], 
	
          RGBColor[1, 0, 0]\  -> \ GrayLevel[0.3], \ 
			
          RGBColor[0. , 0. , 1. ]\  -> GrayLevel[0], 
          RGBColor[0. , 1. , 0.498001]\  -> GrayLevel[0], 
		
          RGBColor[1,  .2,  .4]\  -> GrayLevel[ .9]}, \ 
      DisplayFunction\  -> $DisplayFunction]; \)

... A Queueing Case Study

This case study describes how a real-world situation was turned into a valuable project for students. The students were able to use the WaitingTimes function in a large business meeting to convince the principals that their planned communication links were in error by orders of magnitude.

Revenue Management Services (RMS) is a company specializing in technical assistance, both regulatory and operational, to providers of medical services. One of its ventures involves establishing a wide area computer network (WAN) for processing Medicaid applications for current and prospective nursing home residents. The applications themselves are two or three pages of responses to questions pertaining to a resident’s health and financial situation. However, supporting documentation, including recent bank statements, involves 150 pages on average. The applications with supporting documentation are copied at the nursing home and forwarded to New York City Human Resources Administration (NYC HRA) for a determination of eligibility. For a number of reasons having to do with the sheer volume of the paperwork, the eligibility determinations take months to complete, delaying placement decisions and tying up cash.

RMS has contracted with a software company specializing in optical imaging to develop software that would take a resident’s Medicaid application electronically and attach scanned images of the supporting documentation to a record in a database residing on the network’s server to be located north of NYC. Each nursing home subscribing to this service would have a personal computer (PC) with an attached scanner. The nursing home computers would be attached to the server via one of several possible telecommunication options. Matt, the president of the software company, has recommended a Frame Relay solution that would involve leasing a T1 trunk composed of 24 lines at $980 per month, a $1400 per month access charge, and $40 per month for each nursing home for unlimited calling. With the Frame Relay solution, the nursing home’s computer would never encounter a busy signal, because the communication carrier’s communications computers would merely store information sent by the nursing home until the network server was ready to accept it. This solution is also endorsed by David, RMS’s programmer assigned to this project, and RMS’s marketing managers, Marty and Lew, who like the idea that the customer would never encounter a busy signal.

Tom, a communications consultant, and Karl, a statistical consultant, oppose the Frame Relay solution arguing that is unnecessarily costly for the amount of traffic that the network is likely to see during its first year of operations. Marty has estimated that by the end of the year there will be 40 nursing homes subscribing. Marty and Lew anticipate that the larger homes would tend to be the first to subscribe. The solution favored by Tom and Karl involves using high-speed modems capable of receiving 10 scanned pages per minute at a cost to the nursing home of 8 cents per minute and a cost to RMS of $980 per month.

Compare the two solutions with respect to cost and various waiting line characteristics, such as probability of busy, average number of customers in the system, and so on. State your assumptions and vary the parameters of your model to show under what extreme conditions your recommendation still holds; that is, conduct a sensitivity analysis. Simulate your model using WaitingTimes, which will provide a graphical description of the way the proposed system would operate. Your report should include an executive summary accessible to a nontechnical reader, which explains your figures, graphs, and conclusions.

.  Summary

There are a number of very good software tools available to statisticians for performing data analysis. But Mathematica’s tremendous breadth will generally allow one to do much more. The functions in the VisualStats package illustrated here were quite straightforward to program; readers with better visualization ideas in this area should have no difficulty implementing them in Mathematica. Indeed, statistics is ripe for this sort of work because the community of statisticians has done a lot of excellent work on visualizing data sets. Another area that we have not touched on is the use of the front end to create flashy notebooks for student use. Again, there are a lot of materials on the use of computing in statistics education, and we imagine that many of the activities would be well suited to Mathematica demonstrations.

References

[Feller 1968] Feller, W., An Introduction to Probability Theory and Its Applications, Wiley, New York, 1968.

[Scheaffer 1996] Scheaffer, R. L., Gnanadesikan, M., Watkins, A., and Witmer, J. A., Activity-Based Statistics: Instructor Resources, Springer, New York. 

[Wagon 1995] Wagon, S., Bézier curves, Mathematica in Education and Research, 4:1 (48–53).

About the Authors

Karl Heiner is an associate professor in the Department of Business
Administration at the State University of New York at New Paltz. His doctorate
is from Columbia University, in teaching of college statistics. He uses
Mathematica in teaching introductory statistics courses, in his consulting work,
which includes preparing evidence as an expert witness, and in his research into
more effective estimation in statistical auditing and in the monitoring health
care quality.

Karl Heiner
Department of Business Administration
SUNY at New Paltz
New Paltz, NY12561
kwh@world.com

Stan Wagon uses Mathematica extensively in his research, teaching, and exposition, and is especially appreciative of the new ways Mathematica allows us to look at mathematical objects. He is author or coauthor of several books, including Mathematica in Action, Animating Calculus, VisualDSolve, and Which Way Did the Bicycle Go?. He also teaches a high-altitude course, Rocky Mountain Mathematica, every summer in the mountains of Colorado.

Stan Wagon
Department of Mathematics and Computer Science
Macalester College
St. Paul, MN 55105
wagon@macalester.edu