Sampling Distributions for Random Data

Mathematica in Education and Research
Vol. 4 No.3
Copyright 1995 TELOS/Springer-Verlag

by
Peter Mitic and Peter G. Thomas
Faculty of Mathematics and Computing
The Open University
Milton Keynes, MK7 6AA, UK

Abstract

Mathematica provides the means to study sampling distributions by constructing the necessary theory in a way that emphasises a practical development.  The aim of this paper is to show how Mathematica can be used to study the properties of sampling distributions.  We show how transformations of random numbers can be used to simulate a number of different probability distributions, and illustrate how they may be analysed graphically.   Methods are built up gradually in order to produce increasingly flexible procedures.  A cautionary Mathematica proof of the Central Limit Theorem is also given.

Introduction

The Central Limit Theorem is commonly taught as part of the 16-18 school mathematics curriculum.  The impact of the theorem is considerable, since if the  sample size is sufficiently large, observations from any background population can be treated in exactly the same way.  The impact of the theorem can be a surprise for students who are used to dealing with many types of distribution.  Alternatively, if they have not, the scope of the theorem can be missed.

The Central Limit theorem is often taught in the context of inference, so that the two become inextricably (and wrongly) associated.  The computational difficulties in providing demonstrations of the theorem also cause problems.  Such demonstrations are frequently omitted, either because they are not available or because they are not considered necessary.  We argue that practical demonstrations of the Central Limit Theorem are necessary for any study of this topic in order to acquaint the student with the conceptual framework surrounding the theorem and to engender an appreciation of the context in which it occurs.   Here, we present a discussion of how Mathematica may be used to present the theorem in a particularly simple and efficient way.

We stress that demonstrations of a result, however numerous, cannot constitute a proof. Consequently, the demonstrations discussed here should not replace a proof.  If the mathematics necessary to provide the proof are beyond the scope of the course being taught, the best that can be done is to state the theorem without the proof, noting that the proof has been omitted.  We have, however, developed a Mathematica proof of the theorem, and discuss the advantages and disadvantages of using it.

Qualitative, quantitative and visual approaches

Statements of the theorem, such as given in [DHJ95], are usually in two main parts.

The first states that if the mean and standard deviation of a random variable are m\greek and s\greek respectively, then the corresponding results for the mean of samples of size n are m\greek and (s\greek)/Sqrt[n] respectively.  It is not hard, although rather tedious, to demonstrate these results using a calculator. They are also not difficult to prove.  It is not practical to use a calculator to do either a large number of trials or use a large sample size.  The computer is the practical way out.  Thus, both qualitative and quantitative results can be achieved for the mean and standard deviation of the sampling distribution.

Having generated results, they can then be hard to visualise without extensive sampling, which is also best done by computer.  It is necessary to plot a histogram of the original data and a histogram of means of samples of given size.  The third part of the statement of the Central Limit Theorem, which is that means of samples of given size are Normally distributed, then becomes meaningful.  We give, below, some simple techniques to produce histograms.  Visualisation of the data in this way directly illustrates the impact and scope of the Central Limit Theorem.

The technique of inverting the distribution function for a non-uniform distribution to obtain random numbers which have that non-uniform distribution is used here to provide demonstrations involving non-uniform distributions.  The reader is referred to, for example, [Mitic 1982] for details and discussion of the method.   We will give an example for the normal distribution.  The computer is ideal for applying a transformation to every element in a large set.

In order to stress the programming and pedagogic aspects of this work, we have not used any predefined package functions.  Suitable packages, InverseStatisticalFunctions, for example, are available, and are discussed in [ST94].

Generation and frequency grouping of random data

The primary object with which we shall work is the random number.  This is convenient because Mathematica, in common with most mathematical software, can generate pseudo-random numbers which are theoretically uniformly distributed in the open interval (0,1).  If X is a random variable which has this distribution, we write X ~ U(0,1).  Random numbers also provide a way of visualising a theoretical probability  distribution.

In the following example we investigate a U(0,1) background population.  The following line generates 10 random numbers from this distribution.

data = Table[Random[],{10}]

Having generated some random numbers, we produce a frequency distribution by categorising them into `buckets'.  For example, given a list of random numbers in the range (0,1), we might want to find how many of them are in each of the intervals 0-0.1, 0.1-0.2, ... , 0.8-0.9 and 0.9-1.0.  In order to do this, three procedures are needed.  The first, LeftBoundaries, determines left boundaries for the class intervals, all of which are assumed to have the same width.  The second,  InitialiseCount, is a frequency count, and is initialised to zero for each interval.  Third, the procedure IncrementFrequencies determines the correct `bucket' to increment, given a target input.

LeftBoundaries[xmin_, xmax_, xstep_]:=
Module[{r},
   Table[xmin + r xstep,
        {r,0,(xmax-xmin)/xstep - 1}]]

InitialiseCount[boundaries_]:= 
   Map[{#,0}&, boundaries]

IncrementFrequencies[target_, fc_]:=
Module[{p, lefts, x},
   lefts = Map[First[#]&, fc];
   p = Length[DeleteCases[-lefts+target,
                         x_?Negative]];
   fc /. fc[[p]] :> {fc[[p,1]],fc[[p,2]]+1}]

This example partitions the interval (0,1) into ten equal intervals, as above, generates 100 U(0,1) random numbers and initialises a count.  A Fold construct is used to perform the increment repeatedly for each random number in the list, and illustrates how Fold can be used very efficiently.

data = Table[Random[],{100}];

lb = LeftBoundaries[0,1,.1]

initf = InitialiseCount[lb]

fcount = Fold[IncrementFrequencies[#2,#1]&, 
             initf, data]

Drawing a histogram can be done by replacing each pair in the list fcount by a graphic object which represents a rectangle of width xstep, with base on the horizontal axis.  It is also useful to be able to colour the bars of the histogram so that two histograms can be visibly superimposed at a later stage.

Histogram[fc_, xstep_,{r_,g_,b_}]:=
   Graphics[fc /. {p_,q_?NumberQ} :> 
      {RGBColor[r,g,b],
       Rectangle[{p,0},{p+xstep,q}]}]

Show[Histogram[fcount,.1,{1,0,0}],
    Axes->True,
    AxesLabel->{FontForm["Random number",{"Times-Italic",10}],
                FontForm["Frequency",{"Times-Italic",10}]}]

The Central Limit Theorem for Uniformly Distributed Random Numbers

Having considered a background population, we now turn to sample distributions.

In order to illustrate the Central Limit Theorem, we need to generate a list of means of random numbers.  The function Means, below, generates a table of n means of samples of size ssize.  A histogram of the distribution of these means can then be produced.

Means[data_, ssize_]:= 
Module[{pts},
   pts = Partition[data,ssize];
   Map[((Plus @@ #)/ssize)&,pts]]

data = Table[Random[],{400}];

meandata = Means[data,4];

lb = LeftBoundaries[0,1,.1];

initmeanf = InitialiseCount[lb];

meancount = 
Fold[IncrementFrequencies[#2,#1]&,
    initmeanf, meandata]

Show[Histogram[meancount,.1,{0,0,1}],
    Axes->True,
    AxesLabel->{
    FontForm["Means, sample size 4",{"Times-Italic",10}],
    FontForm["Frequency",{"Times-Italic",10}]}]

This histogram provides the visual overview of the Central Limit Theorem which we seek.  Comparing it with the uniform profile for the background distribution, it is clear that the sampling distribution is very different.  The familiar bell shape of the Normal distribution is apparent, and its form is better defined if a larger sample size is used.  In order to quantify the computation, we can calculate the mean and standard deviation for the sample.  The mean of means corresponds closely with the population mean, 1/2.  Multiplying the sample standard deviation by Sqrt[4] corresponds closely to the theoretical population standard deviation, 1/Sqrt[12].  The theoretical results are obtained from the integrations which follow.  They use the probability density function for the U(0,1) distribution defined by

sampleMean[data_]:= (Plus @@ data)/Length[data]

sampleVariance[data_]:=
   (Plus @@ (data^2))/Length[data] -
   sampleMean[data]^2

sampleSD[data_]:= Sqrt[sampleVariance[data]]

sampleMean[meandata]

sampleVariance[meandata]

sampleSD[meandata]

theoreticalMean =
   Integrate[x 1, {x,0,1}]

theoreticalSD =
   Sqrt[Integrate[x^2 1, {x,0,1}] -
   theoreticalMean^2]

theoreticalSD/Sqrt[4] //N

The numerical values for the sample mean and population mean correspond reasonably well.  So do the numerical values for the sample standard deviation and the population standard deviation divided by Sqrt[4].  This is despite the small sample size.  Several trials like this indicate that there is usually a close correspondence between the two means concerned, but less so with the two standard deviations, unless the sample size is larger (typically about 10).

It is instructive to superimpose the profiles of the background and the sampling distributions.  In all cases, the background profile appears in red and the sampling profile appears in blue (see [ST94]).

data = Table[Random[],{400}];

meandata = Means[data,4];

lb = LeftBoundaries[0,1,.1];

initmeanf = InitialiseCount[lb];

fcount = 
   Fold[IncrementFrequencies[#2,#1]&,
       initmeanf,data]

hb = Show[Histogram[fcount,.1,{1,0,0}],
         Axes->True,
         AxesLabel->{
           FontForm["U(0,1), sample size 4",{"Times-Italic",10}],
           FontForm["Frequency",{"Times-Italic",10}]},
         PlotRange->All,DisplayFunction->Identity]

meancount = 
   Fold[IncrementFrequencies[#2,#1]&,
       initmeanf, meandata]

hs = Show[Histogram[meancount,.1,{0,0,1}],
         Axes->True,
         DisplayFunction->Identity]

Show[{hb,hs},
    DisplayFunction->$DisplayFunction]

Transforming data to another distribution

So far we have used data from a U(0,1) distribution.  If we require data from some other continuous distribution, we can obtain them by transforming the data.  If this cannot be done, an alternative technique must be used.  Such alternatives tend to involve other ways of processing the random numbers generated, and this applies in particular to discrete random variables.  A limited discussion is given below.

For a first example of the transformation technique, we transform to an exponential distribution.  The random variable, X ~ Exponential(mu), has probability density function
$$
some math here
$$

The inverse of its distribution function is Log[1-x]/(-mu). We can implement this transformation as follows.

SetAttributes[ToExp, Listable];

ToExp[x_/;(x>=0 && x<1), mu_]:=
   Log[1-x]/(-mu)

Thus, in order to generate random numbers from an Exponential(2) distribution we can Map the function ToExp to data from a U(0,1) distribution.

U01data = Table[Random[], {5}]

Map[ToExp[#,2]&,U01data]

The following statements for the Exponential(2) distribution parallel the steps required to produce histograms for the U(0,1) distribution and its sampling distribution.

lb = LeftBoundaries[0,1.6,.1];

initexpf = InitialiseCount[lb];

U01data = Table[Random[],{400}];

expdata = Map[ToExp[#,2]&,U01data];

expcount = 
   Fold[IncrementFrequencies[#2,#1]&,
       initexpf,expdata]

hb = Show[Histogram[expcount,.1,{1,0,0}],
         Axes->True,
         AxesLabel->{
         FontForm["Exp(2), sample size 4",{"Times-Italic",10}],
         FontForm["Frequency",{"Times-Italic",10}]},
         PlotRange->All,DisplayFunction->Identity]

expMeandata = Means[expdata,4];

expMeancount = 
   Fold[IncrementFrequencies[#2,#1]&,
       initexpf, expMeandata]

hs = Show[Histogram[expMeancount,.1,{0,0,1}],
         Axes->True,
         DisplayFunction->Identity]

Show[{hb,hs},
    DisplayFunction->$DisplayFunction]

It is clear that this profile has elements of a bell shape, it is skewed to the left.  This reflects the skewness of the background distribution and the small sample size.  The profile looks much more ``normal'' with  a sample size of 10 or more.  The reader is invited to investigate the effect of different sample sizes on the profile of the sampling distribution, and to comment on the question ``how small can a large sample be into order to produce a reasonable normal profile in this context?''

In order to be flexible in adapting our method to other distributions, it is convenient to write a Module, BSGraph, which combines all the necessary operations. Here it is for a symmetric V-shaped distribution with probability density function
$$
some math here
$$

for which the inverse distribution function is
$$
more math here
$$

BSGraph[xmin_,xmax_,w_,invF_,n_,ssize_,xlabel_]:=
Module[{lb,initf,U01data,Tdata,Tcount,hb,hs,
        Meandata,MeanCount},
   lb = LeftBoundaries[xmin,xmax,w];
   initf = InitialiseCount[lb];
   U01data = Table[Random[],{n}];
   Tdata = Map[Nest[invF,#,1]&,U01data];
   Tcount = Fold[IncrementFrequencies[#2,#1]&, 
                initf,Tdata];
   hb = Show[Histogram[Tcount,w,{1,0,0}],
            Axes->True,
            AxesLabel->{xlabel,"Frequency"},
            PlotRange->All,DisplayFunction->Identity];

   Meandata = Means[Tdata,ssize];
   Meancount = Fold[IncrementFrequencies[#2,#1]&,
                   initf, Meandata];
   hs = Show[Histogram[Meancount,w,{0,0,1}],
            Axes->True,
            DisplayFunction->Identity];
   Show[{hb,hs},DisplayFunction->$DisplayFunction]]

invF[x_ /; 0<=x<=1/2]:= (1-Sqrt[1-2x])/2

invF[x_ /; 1/2<x<=1]:= (1+Sqrt[2x-1])/2

invF[x_]:=0 (* elsewhere *)

BSGraph[0,1,.1,invF[#]&,400,4,"V, sample size 4"]

Mathematica is particularly well suited to deal with a piecewise continuous function like this because of the natural way in which such a function may be defined.  This example is particularly noteworthy for the way in which the highly skewed V-shaped background distribution gives rise to a recognisable normal-shaped sampling distribution, even with the small sample size of 4.  A sample size of 2 or 3 is interesting.

Sampling from a Normal distribution

In order to produce random numbers which are normally distributed, we require a more sophisticated transformation of U(0,1) random numbers.  We have adapted a method, details of which may be found in [HH64].  An alternative can be found in [AS72].  A  normally distributed random variable with mean 0 and variance 1, denoted N(0,1), may be simulated by the function NormalRandom[].  Theoretically, the results can vary from -infinity to infinity. This causes problems if we try to use the Fold function with LeftBoundaries, so we truncate the domain to [-2.5,2.5].  Consequently, this is not a true random sample from a normal distribution.  The probability of obtaining a N(0,1) observation outside this range is small, and we can afford to ignore this region.

NormalRandom[]:= Module[{w,w1,w2,n1,n2}, 
   w = 1;
   While[w >= 1,
        w1 = 2 Random[]-1;
        w2 = 2 Random[]-1;
        w = w1^2 + w2^2
        ];
   n1 = N[w1 Sqrt[-2 Log[w]/w]];
   n2 = N[w2 Sqrt[-2 Log[w]/w]];

(* Either n1 or n2 can be returned.
   Truncate the domain first *)

   n2 = If[n2 > 2.5, 2.5, n2];
   n2 = If[n2 < -2.5, -2.5, n2];
   n2]

It works in the same way as the built-in Random.

Table[NormalRandom[], {4}]

The procedure BSGraph cannot, as presented above, use NormalRandom.  In order to make this possible, we have to provide an option to use an alternative to the default function, Random.  The method of doing this is explained in [Mae91].  Here is a redefined version of BSGraph which includes this feature.

Options[BSGraph] = {RandomNumberGenerator->Random}

BSGraph[xmin_,xmax_,w_,invF_,n_,ssize_,xlabel_,opts___]:=
Module[{lb,initf,U01data,Tdata,Tcount,hb,hs,
        Meandata,MeanCount,RNGen},
   RNGen = RandomNumberGenerator /. {opts} /.
              Options[BSGraph];
   lb = LeftBoundaries[xmin,xmax,w];
   initf = InitialiseCount[lb];
   U01data = Table[RNGen[],{n}];
   Tdata = Map[Nest[invF,#,1]&,U01data];
   Tcount = Fold[IncrementFrequencies[#2,#1]&, 
                initf,Tdata];
   hb = Show[Histogram[Tcount,w,{1,0,0}],
            Axes->True,
            AxesLabel->{xlabel,"Frequency"},
            PlotRange->All,DisplayFunction->Identity];

   Meandata = Means[Tdata,ssize];
   Meancount = Fold[IncrementFrequencies[#2,#1]&, 
                   initf, Meandata];
   hs = Show[Histogram[Meancount,w,{0,0,1}],
            Axes->True,
            DisplayFunction->Identity];
   Show[{hb,hs},
       DisplayFunction->$DisplayFunction]]

Using the RandomNumberGenerator->NormalRandom option produces curious spiky profiles which are unlike the corresponding profiles produced by other software.  This is partly due to the small interval width in this case, and partly due to the random number generator.

invF[x_] := x (* identity function *)

BSGraph[-2.5, 2.5, .1, invF[#]&,
       400,4, "N(0,1), sample size 4",
       RandomNumberGenerator -> NormalRandom]

A Discrete Distribution: the Binomial

Having already catered for alternative random number generators, it is easy to simulate discrete events.  Here is an example of a binomial random variable.  A general binomial function, which takes two parameters, one for the number of independent trials and the other for the probability of success of any given trial, is defined first.  However, a function which takes no arguments is required for BSGraph, so we also have to define a specific binomial function.

BinomialRandom[p_,n_] := Module[{ct,i},      
   ct = 0;
   For[i=1, i<=n, i++,
      If[Random[] < p, ct = ct + 1]
      ];
   ct]

BinRandom[]:= BinomialRandom[.25,21]

A check shows that it behaves reasonably, and the graphs are as expected.

Table[BinRandom[], {15}]

invF[x_]:=x (* identity function *)

BSGraph[0,21,1,invF[#]&,400,4,
   "Bin(.25,21), sample size 4",
   RandomNumberGenerator->BinRandom]

Other distributions

In  [MT95] we discuss software implemented in Visual Basic which uses a number of transformation techniques for random numbers to simulate some common discrete distributions, including the Binomial and Poisson.  In the same article, we also use a Monte Carlo technique to simulate sampling from some more awkward distributions.

We invite the reader to apply the same, and other, techniques using Mathematica.

Fitting a normal curve

It is possible to fit a normal curve to the sampling distributions we have generated.  We first define a normal curve with parameters m and s, as defined by the mean and standard deviation of the data.  This needs to be scaled to account for the total frequency and the interval widths.  In our example, which uses a U(0,1) background distribution, we use 100 data and group them into intervals, each of width 0.1.  Hence, the scale factor is 0.1 times 100 = 10.  This reflects the fact that the total area of the histogram is 10, and not 1.

Phi[x_, m_, s_]:=(E^(-(x-m)^2/(2s^2)))/(s Sqrt[2Pi])

data = Table[Random[],{400}];

Meandata = Means[data,4];

lb = LeftBoundaries[0,1,.1];

initmeanf = InitialiseCount[lb];

Meancount = 
   Fold[IncrementFrequencies[#2,#1]&,
       initmeanf, Meandata]

A histogram is combined with the scaled normal curve.

sm = sampleMean[Meandata]

sv = sampleVariance[Meandata]

hist = Show[Histogram[Meancount,.1,{0,0,1}],
           Axes->True,
           AxesLabel->{
           FontForm["U(0,1), sample size 4",{"Times-Italic",10}],
           FontForm["Frequency density",{"Times-Italic",10}]},
           DisplayFunction->Identity]

nplot = Plot[10 Phi[x,sm, Sqrt[sv]], 
            {x,0,1},
            AxesOrigin->{0,0},
            PlotStyle->{{RGBColor[0,1,0]}},
            DisplayFunction->Identity]

Show[{hist, nplot},
    DisplayFunction->$DisplayFunction]

Proof

Finally, we give a proof of the Central Limit Theorem, with some caveats.  It is an example of a proof which needs a lot of setting up, and once this has been done, is very short with the aid of Mathematica.  The usual proof uses the moment generating function of a random variable.  If a random variable, X, has a probability density function, f(x), which is positive on a subset, S, of the real line, then the moment generating function of X, M(t), can be obtained by integrating exp(x t) f(x) with respect to x on S.

As a preliminary, we establish the result that the moment generating function for a N(0,1) random variable is exp(-(t^2)/2).

MGFNormal[t] =
Integrate[E^(x t) Phi[x,0,1],
         {x, -Infinity, Infinity}]

In order to proceed, we must set up a formal statement of the Central Limit Theorem.  Mendenhall (in [MSW86]) gives a particularly clear statement.

Let Y1, Y2, ..., Yn be independent and identically distributed random variables, each
with mean m and variance s^2.  Then Ybar = (Y1+Y2+ ... + Yn)/n  is approximately Normally distributed with mean m\greek and variance ((s\greek)^2)/n for large n.

Standardising each random variable Yi by defining Zi = (Yi - m)/s,  we have that the mean of Zi is 0 and the variance of each Zi is 1.   The moment generating function of each Zi can then be written as M(t) = 1 + (t^2)/2 + t^3 R(t)  where R(t) is a power series in t which comprises moments of Zi  of degree 3 or more.  This is the starting point for the Mathematica session.  We then define Un as below and express it in terms of the Zi.
$$
math here
$$

Since the Yi are independent, so are the Zi.  Hence, the moment generating function of Un is the product of the moment generating functions of the Zi/Sqrt[n].  Taking a limit as n tends to infinity gives the moment generating function of a Normal(0,1) random variable.  Given that the moment generating function for any given random variable is unique, this shows that Un has a Normal(0, 1) distribution.  Hence, the non-standardised random variable, Ybar, has a Normal(m, (s^2)/n) distribution.

The result should therefore arrive in two lines of input, but, Mathematica cannot obtain the limit for the general remainder function, R(t) unless we specify that R(t) is a polynomial in t.   We really need this important stage, because it takes rather more analysis to do it by hand.  Here is how it may be done.

M[t_]:= 1 + t^2/2 + t^3 R[t]

Limit[M[t/Sqrt[n]]^n,n->Infinity]

R1 /: PolynomialQ[R1,t]=True

Mpoly[t_]:= M[t] /. R[t]-> R1

Limit[Mpoly[t/Sqrt[n]]^n ,n->Infinity]

Using this method, we have achieved in three lines what Mendenhall [MSW86] takes most of a page to do.

Finally, here is the caveat.  It is easy to produce a 'proof' which gives the correct result from a wrong premise.  We discuss this, and related problems, in more detail in [MT94].  If the approximation for R(t) is too harsh, the same result is nevertheless obtained.  Two examples are R(t) = a, a constant, and R(t) = 0.

Mconst[t_]:= M[t] /. R[t]-> a

Limit[Mconst[t/Sqrt[n]]^n,n->Infinity]

Mzero[t_]:= M[t] /. R[t]-> 0

Limit[Mzero[t/Sqrt[n]]^n,n->Infinity]

Right, but wrong!

Conclusion

The impact and importance of the Central Limit Theorem can be easily lost without a visual approach, which complements probabilistic and statistical calculations.  Using Mathematica, we have been able to simulate number of different probability distributions, generate samples and analyse them graphically.  Mathematica is undoubtedly slower than a specialist statistics package, but has the advantage of an integrated environment for doing algebraic manipulation in the context of sampling.  Students can also benefit from building the theory from scratch rather than using a ``black box'' approach.

Two qualitative conclusions arise from this work.  First, about 100 observations of the mean are sufficient in most cases.  The shapes of the background and of the sampling distributions are reasonably clear, especially if more buckets are used.  Secondly, a "large" sample size is about 10, provided that the background population is not very skewed.  We used a sample size of 4 here to obtain results faster, and sampling distributions are skewed if the background distribution is skewed.

We stress that the only quantitative results obtained here are trials in which the mean and standard deviation of a background population and a corresponding sampling population are calculated.  An algebraic approach is important, which is why we presented a proof as well.  However, first courses in statistics often omit the proof because it is mathematically demanding.  In such circumstances we recommend that it is made clear that a result has been stated without proof, and that none of the demonstations decribed here can replace a proof.

References

[AS72]	M. Abramowitz and I. Stegun.  Handbook of Mathematical Functions. 1972.  Dover.

[DHJ95] F. Daly, D. J. Hand, M.C. Jones, A.D. Lunn and K.J. McConway. Elements of Statistics. 1995. Addison Wesley.

[HH64] J.M. Hammersley and D.C.Handscomb.  Monte Carlo Methods. 1964. Chapman and Hall.

[Mae91] R.Maeder.  Programming in Mathematica. 1991.  Addison Wesley.

[Mit82] I. Mitrani.  Simulation Techniques for Discrete Event Systems. 1982. Cambridge University Press. 

[MSW86] W. Mendenhall, R. Schaeffer, and D. Wackerly.  Mathematical Statistics with Applications. 1986.  PWS, Boston.

[MT94] P. Mitic and P.G. Thomas.  Pitfalls and limitations of computer
algebra, in Computers and Education 22, 4.  1994	Elsevier Science.

[MT95] P.Mitic and P.G.Thomas. The Central Limit Theorem---Visualised, submitted to Teaching Mathematics and its Applications. 1995.  Institute of Mathematics and its Applications, Southend, UK.	

[ST95]	W. Shaw and J. Tigg.  Applied Mathematica: Getting Started, Getting it Done.  1994.  Addison Wesley.

About the Authors

Peter Mitic and Peter G. Thomas
Faculty of Mathematics and Computing
The Open University
Milton Keynes, MK7 6AA, UK
P.Mitic@soton.ac.uk 
P.Thomas@soton.ac.uk