IMSE Chapter 10

IMSE-10(Statistics).nb
Statistics in Mathematica
RCH 12/31/01
revised 6/26/02 - rjerwin

This notebook explores some of the statistical capabilities of Mathematica.

Descriptive Statistical Operations

In this section we will explore various ways of setting up statistical operations in Mathematica.  

First, we will assume that we have some data stored (as usual in Mathematica) in a list and we want to find various statistical measures of those data such as the mean (average), variance (standard deviation), etc.  Mathematica has convenient built-in functions to evaluate these statistical quantitites.  They are part of the statistical package Statistics`DescriptiveStatistics`, which, of course, must be loaded.

<< Statistics`DescriptiveStatistics`

Now let's make up some data representing, for example the heights (in meters) of a group of people.

heightList = {1.53, 1.77, 1.85, 2.02, 1.94, 2.11, 1.85, 
          1.42};

We can now have Mathematica compute the mean, the variance, and the standard deviation of the data.

meanHeight = Mean[heightList];
heightVariance  =  
        Variance[heightList];
deviation = 
        StandardDeviation[heightList];
Print["meanHeight = ",  meanHeight, ", heightVariance = ", 
      heightVariance,  ", STD = ", deviation]

Mathematica has a nice function called DispersionReport that gives you several measures of the "spread" in your data..

DispersionReport[heightList]

Statistical Properties of Histograms and Other "Frequency" Distributions

In this section, we show how to calculate the mean, standard deviation, etc. for data that represent so-called "frequency" distributions.  For example, we may have data that tell us how many people have various heights (rather than just a list of heights as we had in the example above).  To be specific, let's assume that the data appear in a list with the heights in meters and the corresponding number of people paired off as shown in the next input cell.

dataList = {{1.5, 10}, {1.6, 15}, {1.7, 12}, {1.8, 20}, {1.9, 
            25}, {2.0, 19}, {2.1, 5}, {2.2, 2}};
ListPlot[dataList, AxesLabel -> {"height", "number"}, 
        PlotStyle -> PointSize[ .02], 
        AxesOrigin -> {1.45, 0}];

An aside:  This is the perfect sort of thing to put in that old elementary school standby -- a bar chart!  First we need to load the Graphics package.  Next, the command GeneralizedBarChart requires a third set of numbers in dataList specifying the width of each bar.  Here I've set all of them at 0.05.

<< Graphics`Graphics`
databarList = {{1.5, 10,  .05}, {1.6, 15,  .05}, {1.7, 
            12,  .05}, {1.8, 20,  .05}, {1.9, 25,  .05}, {2.0, 
            19,  .05}, {2.1, 5,  .05}, {2.2, 
            2,  .05}};
GeneralizedBarChart[databarList, AxesOrigin -> {1.45, 0}, 
        AxesLabel -> {"height", "number"}];

Now suppose that we want to find the mean and standard deviation of this distribution.  To calculate the mean, we multiply each height value by the number of people with that height, sum the results, and then divide by the total number of people.

To carry out these operations, we first need to separate the lists.  This is done with the Part[] function in Mathematica.  See IMSE-05(Lists-Arrays) and section 1.8.4 of the Mathematica Book.

\(heightDataList = 
      Part[dataList, All, 
        1] (*Pulls\ out\ the\ first\ "column"\ of\ \(\(data\)\(.\)\)*) \)
\(numberDataList = 
      Part[dataList, All, 
        2] (*Pulls\ out\ the\ second\ "column"\ of\ \(\(data\)\(.\)\)*) \)

Next we calculate the mean and standard deviation.  The operation Apply[Plus,list] (alternatively, you could say Plus@@list) gives the sum of the elements in a list.

totalNumber = 
        Apply[Plus, numberDataList];
meanHeight = 
        Apply[Plus, numberDataList*heightDataList]/
          totalNumber;
Print["totalNumber = ", totalNumber,  ", meanHeight = ", 
      meanHeight]

You should find that the mean height is 1.81 meters (to two significant figures) for this group of 108 people.

Next we compute the standard deviation from the mean.  Recall that the standard deviation is the square root of the average squared-deviation from the mean.  So, first we calculate a list containing the difference of the height from the mean height.  We then square those differences, find the average of those quantities, and finally take the square root.  (We won't fuss over the difference between the standard deviation of a sample and the standard deviation of a population from which the sample is drawn.) 

 If you want to see what the intermediate lists look like, remove the semi-colons at the ends of the input statements.

heightDiffList = 
        heightDataList - meanHeight;
heightDiffSqList = 
        heightDiffList*heightDiffList;
prodDiffList = 
        heightDiffSqList*numberDataList;
averageSqDev = 
        Apply[Plus, prodDiffList]/totalNumber;
stndDevMean = (Sqrt[averageSqDev]);
Print["variance = ", averageSqDev,  ", standard dev. = ", 
      stndDevMean]

So, we see that for this height distribution the mean height is 1.81 meters and the standard deviation from the mean is 0.18 meters (both given to two signficant figures.)  If the height distribution were described by a Gaussian ("normal") distribution, then we would expect to find about 68% of the heights to be within one standard deviation on either side of the mean.  In this example, 84% of the heights lie between 1.6 m and 2.0 m.  That result tells us that our distribution is not a Gaussian distribution, which should have been obvious from the graph given above.

Using Mathematica's Probability Distributions

Mathematica has many handy built-in probability distribution functions.  To illustrate some of these features, let's construct a Gaussian distribution and then analyze it using methods developed in the previous section.  These distribution functions (and related operations) are found in the package Statistics`ContinuousDistributions`.

<< Statistics`ContinuousDistributions`

This Gaussian distribution will have a mean of 1.81 and standard deviation of 0.18 to match what we found above.  The operation PDF[dist,x] gives the numerical value of the distribution at the point x.

ourDist = 
        NormalDistribution[1.81, 
          0.18];

  
xList = Table[x, {x, 0, 3.6, 0.01}];
pdList = 
        Table[PDF[ourDist, x], {x, 0, 3.6, 
            0.01}];
ListPlot[Transpose[{xList, pdList}], PlotJoined -> True, 
        PlotRange -> All, 
        AxesLabel -> {"x", "probability density"}];

Now let's use the calculations developed to calculate the mean and standard deviation.  First, the mean.

total = Apply[Plus, pdList];
meanx = Apply[Plus, pdList*xList]/total;
Print["total = ", total,  ", mean of x = ", meanx]

Now, the standard deviation.

xDiffList = xList - meanx;
xDiffSqList = xDiffList*xDiffList;
prodxList = xDiffSqList*pdList;
varx = 
      Apply[Plus, prodxList]/
        total;   
stndDevx = (Sqrt[varx]);
Print["variance = ", varx,  ", standard deviation = ", 
      stndDevx]

We see that our results agree very closely with the mean and standard deviation of the original distribution.

An aside:  If you want to see the analytic form of the probability distribution function, you can define the distribution with symbols.

PDF[NormalDistribution[μ, σ], x]