MA206 Supplemental Material

Notes in Preparation of a Reading on Empirical Distribution Functions
2004-00-03

LTC Mark Wroth

Introduction and Purpose

This Mathematica notebook is a collection in process notes and calculations related to creating supplemental reading on empirical distribution functions for MA206.  It is used to record notes on the subject in the process of creating the reading material, and to create figures for inclusion in the reading if the Mathematica code is not germain to the reading.

Data Manipulation

To extract a column of data from an array (e.g. one read in as paired data), use the Column or ColumnTake functions (both of which are part of the Statistics`DataManipulation` package).

Ordinary operations can be applied directly to lists.  List sums can be computed with the Sum command.  So to take the sum of the ((Subscript[x, i])^2)from the set read in in datain above

datain

datain[[1, 2]]

xi  =  Column[datain, 1]

Sum[xi[[i]], {i, Length[xi]}]

These functions are particularly useful in computing regression values not output by the standard package.

Empirical Distribution Functions

Overview

This section discusses the mechanics of plotting and manipulating Empirical Distribution Functions (EDF) in Mathematica.

Since this subject is best discussed with the aid of an example, we start by creating the array of data for which we want to plot the EDF. This example is taken from Project II as presented during STAP 04. I have significantly reduced the number of replications in order to be able to show the data without taking too much space.

n = 10;

T = 
      RandomArray[ExponentialDistribution[1. /14], n] + 
        RandomArray[GammaDistribution[2, 14], n] + 
        RandomArray[UniformDistribution[6, 12], n]

Now what we need to do is to create a matrix with the sorted values of the sample each appearing twice in the first column, and the frequency data appearing in the second column.  The exact organization of the data is intended to cause the plotting routines to plot it as a step function.

First, create a vector (list) of the frequency values.  This list consists of the values {0, 1/n, 1/n, 2/n, 2/n, 3/n, 3/n, ... , (n-1)/n, (n-1)/n, 1}.  We create it by first creating a list of the values {1/n, 2/n, ..., 1}

Range[1.0, n]/n

Then we double the list and insert the value 0, converting the nested lists into a single flat list

Flatten[{0, Range[1, n]/n, Range[1, n]/n}]

Notice that every value appears twice except 0.  Now we sort the list and remove one of the 1's at the top of the list, assigning the result to a symbol "freq".

freq = 
      Delete[Sort[Flatten[{0, Range[1, n]/n, Range[1, n]/n}]], 2  n + 1] // 
        N

Each value in the sample also appears twice, and needs to be sorted;

Sort[Join[T, T]]

If we simply make a list of the two lists, we will have a matrix with the right values in it.  However, the list will consist of two rows and n columns, while we want it to be two columns and n rows.  The Transpose function takes care of this.

EDF = Transpose[{Sort[Join[T, T]], freq}]

Now we can plot the resulting pairs.  The "PlotJoined->True" option causes Mathematica to join the plotted points with line segments, which is the desired form of the plot.

EDFplot  =  ListPlot[EDF, PlotJoined -> True]

Other Examples

A uniform random variable

n2  =  250
T2 = 
        RandomArray[UniformDistribution[5, 10], 
          n2];
freq2 = 
        Delete[Sort[Flatten[{0, Range[1, n2]/n2, Range[1, n2]/n2}]], 
            2  n2 + 1] // N;
EDF2 = 
        Transpose[{Sort[Join[T2, T2]], freq2}];
EDFplot2  =  ListPlot[EDF2, PlotJoined -> True]

n3  =  250
T3 = 
        RandomArray[NormalDistribution[10, 5], 
          n3];
freq3 = 
        Delete[Sort[Flatten[{0, Range[1, n3]/n3, Range[1, n3]/n3}]], 
            2  n3 + 1] // N;
EDF3 = 
        Transpose[{Sort[Join[T3, T3]], freq3}];
EDFplot3  =  ListPlot[EDF3, PlotJoined -> True]

n4  =  250
T4 = 
        RandomArray[ExponentialDistribution[0.5], 
          n4];
freq4 = 
        Delete[Sort[Flatten[{0, Range[1, n4]/n4, Range[1, n4]/n4}]], 
            2  n4 + 1] // N;
EDF4 = 
        Transpose[{Sort[Join[T4, T4]], freq4}];
EDFplot4  =  ListPlot[EDF4, PlotJoined -> True]

Fitting A Model to the EDF

Once the EDF has been constructed, we need to fit a model of some sort to it.  The fundamental process is to select parameters for the model that cause it to fit "well" to the data; there are several ways to approach this problem.

Estimator Approach

One approach to this problem is one that fits (pun intended) well with traditional statistics: define an estimator for the parameter(s) in question and use the point estimates resulting from it.  In this case, we will try to fit the Gamma distribution, finding point estimates for E(X) = αβ and V(X) = α(β^2)

LocationReport[Column[EDF, 1]]

Variance[Column[EDF, 1]]

Solve[{α β  ==  
          Mean[Column[EDF, 1]], α  (β^2)  ==  
          Variance[Column[EDF, 1]]}, {α,  β}]

These are our estimates for α and β; check to make sure they really have the desired properties (which is mostly a check that we did the solution right).

Mean[
      GammaDistribution[7.770079564724459`, 6.468184202701226`]]

Variance[
      GammaDistribution[7.770079564724459`, 6.468184202701226`]]

Plot the fitted graph

ListPlot[
      Table[{Quantile[
            GammaDistribution[7.770079564724459`, 6.468184202701226`], i], 
          i}, {i, 0,  .999,  .001}]]

And now the fit and the EDF on the same axis.

MultipleListPlot[
      Table[{Quantile[
            GammaDistribution[7.770079564724459`, 6.468184202701226`], i], 
          i}, {i, 0,  .9999,  .02}], EDF, PlotJoined -> True]

Clearly this is a reasonably good fit of the distribution; since it's part of Project II, let's apply the same technique to an exponential.

MultipleListPlot[
      Table[{Quantile[
            GammaDistribution[7.770079564724459`, 6.468184202701226`], i], 
          i}, {i, 0,  .9999,  .02}], 
      Table[{Quantile[ExponentialDistribution[1/Mean[Column[EDF, 1]]], i], 
          i}, {i, 0,  .9999,  .02}], 
EDF, 
      PlotJoined -> True]

Clearly, while it may have the same mean, the exponential distribution is not as good a model of the simulated data as the Gamma distribution.

Least Squares Fit

Quantile[GammaDistribution[2, 14], Column[EDF, 2]]

((Column[EDF, 1] - 
              Quantile[GammaDistribution[2, 14], 
                Column[EDF, 2]])^2)[[1]]

Compute the sum of squared errors.  Note that in this case we are deleting the last value from the sum, as the Gamma is not well defined at that point.

Sum[((Column[EDF, 1] - 
                Quantile[ExponentialDistribution[1. /14], 
                  Column[EDF, 2]])^2)[[i]], {i, 1, 
        Length[Column[EDF, 1]] - 1}]

NMinimize[{Sum[((Column[EDF, 1] - 
                    Quantile[GammaDistribution[α, β], 
                      Column[EDF, 2]])^2)[[i]], {i, 1, 
            Length[Column[EDF, 1]] - 1}]}, {α, β}]

fitted  =  
      ListPlot[Transpose[{Quantile[ExponentialDistribution[λ], 
                Delete[Column[EDF, 2], 2  n]],  
              Delete[Column[EDF, 2], 2  n]}]] /. λ -> 
          3.8412492407974623`

2  n

Transpose[{Quantile[ExponentialDistribution[λ], 
            Delete[Column[EDF, 2], 2  n]],  
          Delete[Column[EDF, 2], 2 *n]}] /. λ -> 
        3.8412492407974623`

ListPlot[{{0.`, 0.`}, {0.02742871109190325`, 
          0.1`}, {0.02742871109190325`, 0.1`}, {0.05809140134522591`, 
          0.2`}, {0.05809140134522591`, 0.2`}, {0.09285389246564094`, 
          0.3`}, {0.09285389246564094`, 0.3`}, {0.132984243339529`, 
          0.4`}, {0.132984243339529`, 0.4`}, {0.1804483742419288`, 
          0.5`}, {0.1804483742419288`, 0.5`}, {0.2385397755871547`, 
          0.6`}, {0.2385397755871547`, 0.6`}, {0.31343261758145774`, 
          0.7`}, {0.31343261758145774`, 0.7`}, {0.41898814982908356`, 
          0.8`}, {0.41898814982908356`, 0.8`}, {0.5994365240710123`, 
          0.9`}, {0.5994365240710123`, 0.9`}}]

Percentiles

The fundamental definition of a percentile is the value of a random variable such that the defined percent of the distribution lies below the stated value. Devore's notation is not necessarily helpful in understanding this concept.

In order to find a percentile, a direct resort to first principles can be effective.  For the Standard Normal distribution, the PDF is

PDF[NormalDistribution[0, 1], x]

In order to place, say, 0.9 below the percentile c, this equation must be true.

(∫_(-∞)^c) ((E^(-((x^2)/2)))/(Sqrt[2 π])) \[DifferentialD]x == 0.9

If we could solve this equation for c, it would be the desired percentile.  However, this equation is not trivial to solve ... by hand.

Solve[(∫_(-∞)^c) ((E^(-((x^2)/2)))/(Sqrt[2 π])) \[DifferentialD]x == 0.9, c]

The Solve function is not always the best approach to this, however.  Depending on the form of the PDF, Solve may be unable to find a solution (e.g. with the (χ^2)distribution).

PDF[ChiSquareDistribution[ν], x]

Solve[(∫_(-0)^c)(( (2^(-ν/2)) (E^(-x/2)) (x^(-1 + (ν/2))))/Gamma)[(ν/2)] \[DifferentialD]x == 0.9, c]

Getting rid of the Γ function by evaluating it (possible in the specific case, if not the general case) does not appear to solve the problem.

Solve[(∫_0^c)(( (E^(-x/2)) (x^(13/2)))/(135135 (Sqrt[2 π]))) \[DifferentialD]x == 0.9, 
        c] // N

An alternative approach (suggested by Andy Glenn) is to use the FindRoot function:

?FindRoot
"FindRoot[lhs==rhs, {x, x0}] searches for a numerical solution to the equation lhs==rhs, starting with x=x0. FindRoot[{eqn1, eqn2, ... }, {{x, x0}, {y, y0}, ... }] searches for a numerical solution to the simultaneous equations eqni." More…

FindRoot[(∫_0^c)(( (2^(-ν/2)) (E^(-x/2)) (x^(-1 + (ν/2))))/Gamma)[(ν/2)] \[DifferentialD]x == 0.9, {c, 10}] /. ν -> 15

There is an easier way, at least in that Mathematica defines the appropriate operation:

Quantile[NormalDistribution[0, 1],  .9]

Quantile[ChiSquareDistribution[15], 0.9]

Authoring Tools

Tools to help with maintaining documents

<< AuthorTools`

Notebooks[]

nb = %[[1]]

NotebookName[nb]

References And Further Reading

Abell, Martha L..  Statistics with Mathematica.  San Diego, CA: Academic Press, 1998.  A reference book on statistical functions and operations in Mathematica, using both pre-defined and author-defined packages.  Includes a CDROM of related material.

Devore, Jay L..  Probability and Statistics for Engineering and the Sciences, 6th ed.. Belmont, CA: Thomson (Brooks/Cole), 2004.  The course text for MA206 in 2004; the Department has been using previous editions of this text for a number of years.

Gray, John W.  MASTERING MATHEMATICA: Programming Methods and Applications, 2nd Ed. 
ACADEMIC - ELSEVIER, 1998. Contents: Quick Trip Through Elementary Mathematics; Interacting with Mathematica; More About Numbers & Equations; Built-In Graphics & Sound; Mathematica Language; Functional Programming; Rule Based Programming; Procedural Programming; Object-Oriented Programming; Graphics Programming; Polya Pattern Analysis; Object-Oriented Graph Theory; Differentiable Mappings; Critical Points & Minimal Surfaces. CD-ROM INCLUDED.
ACADEMIC - ELSEVIER S ISBN: 0122961056   PGS: 629 List: $49.95     YOUR PRICE: $47.45 From the OpAmp web site and  potentially of interest to me as a student of both mathematics and a teacher of statistics.  Descriptions, prices, and availability are as of 2004-04-11.

Heiner, Karl and Stan Wagon."Statistics in the Classroom and the Courtroom".MATHEMATICA IN ACTION for Issue 8-1 of Mathematica in Education and Research. (http://library.wolfram.com/infocenter/Articles/1730/)

Wagon,Stan, Mathematica in Action,2/e,2nd Ed. SPRINGER VERLAG, 1999.
"Mathematica in Action, 2nd Edition," is designed both as a guide to the extraordinary capabilities of Mathematica as well as a detailed tour of modern mathematics by one of its leading expositors,Stan Wagon.Ideal for teachers,researchers,mathematica enthusiasts.This second edition of the highly sucessful W.H.Freeman version includes an 8 page full color insert and 50% new material all organized around Elementary Topics,Intermediate Applications,and Advanced Projects.In addition,the book uses Mathematica 3.0 throughtout.Mathematica 3.0 notebooks with all the programs and examples discussed in the book are available on the TELOS web site (www.telospub.com).These notebooks contain materials suitable for DOS,Windows,Macintosh and Unix computers.Stan Wagon is well-known in the mathematics (and Mathematica) community as Associate Editor of the "American Mathematical Monthly," a columnist for the "Mathematical Intelligencer" and "Mathematica in Education and Research," author of "The Banach-Tarski Paradox" and and "Unsolved Problems in Elementary Geometry and Number Theory (with Victor Klee), as well as winner of the 1987 Lester R. Ford Award for Expository Writing.
SPRINGER VERLAG 	H 	ISBN: 0387982523 	  PGS: 608 	List: $77.95     YOUR PRICE: $74.05 From the OpAmp web site and  potentially of interest to me as a student of both mathematics and a teacher of statistics.  Descriptions, prices, and availability are as of 2004-04-11.

Stephen Wolfram.The Mathematica Book, 5th ed.. Champaign, IL: Wolfram Media/Cambridge University Press, 2003.  This text on Mathematica is available both in hard copy and (perhaps more importantly) as part of the Mathematica on-line help system.

Mathematica Definitions

The following definitions and package commands have accumulated as needed for various lessons. In general I have added packages and definitions as needed, without removing any.  So this set may include packages or commands that are not needed in today's lesson.  However, at least on the computers I've been working on, having the additional commands has not been a problem.

Note that these cells are normally marked as "initialization cells".  This causes Mathematica to try to ensure they are evaluated before the main cells in the worksheet are evaluated.

<< Statistics`ContinuousDistributions`

<< Statistics`DescriptiveStatistics`

χcrit[α_, ν_]  :=  
      Quantile[ChiSquareDistribution[ν], 1 - α]

<< Statistics`ConfidenceIntervals`

<< Graphics`FilledPlot`

<< Graphics`MultipleListPlot`

<< Statistics`DataManipulation`

showPDF[dist_, ll_, ul_, lpl_, upl_]  :=  
        Show[{Plot[PDF[dist, x], {x, lpl, upl}, 
              DisplayFunction -> Identity], 
            FilledPlot[PDF[dist, x], {x, ll, ul}, 
              DisplayFunction -> Identity]}, 
          DisplayFunction -> $DisplayFunction];

showPDF[NormalDistribution[0, 1], -1, 2, -4, 4]