Introduction to Statistical Computing with Mathematica

Darren Glosemeyer

Wolfram Research, Inc.

Brief Introduction

Built-in Functions

There are a number of functions built into the Mathematica 5.1 kernel of particular relevance to statistics.  These include descriptive statistics functions such as Mean, Variance, StandardDeviation, Median, and Quantile, and fitting functions such as Fit and FindFit.

There are also a number of standard add-on packages included specifically for statistical computations, including descriptive statistics, univariate and multivariate distributions, linear and nonlinear regression, statistical plots, and cluster analysis.

({{"Statistics`ANOVA`"}, {"Statistics`ClusterAnalysis`"}, {"Statistics`ConfidenceIntervals`"}, {"Statistics`ContinuousDistributions`"}, {"Statistics`DataManipulation`"}, {"Statistics`DataSmoothing`"}, {"Statistics`DescriptiveStatistics`"}, {"Statistics`DiscreteDistributions`"}, {"Statistics`HypothesisTests`"}, {"Statistics`LinearRegression`"}, {"Statistics`MultiDescriptiveStatistics`"}, {"Statistics`MultiDiscreteDistributions`"}, {"Statistics`MultinormalDistribution`"}, {"Statistics`NonlinearFit`"}, {"Statistics`NormalDistribution`"}, {"Statistics`StatisticsPlots`"}})

Of course, there are lots of other symbolic, numeric and graphic functions in Mathematica that are useful for statistical computations, though not specific to statistics.

A Mathematical Statistics Example

<< Statistics`

As an example of symbolic statistical computation, let's consider the first four moments of the (χ^2) distribution.  These can be computed directly using the ExpectedValue function.

evs = 
      Table[ExpectedValue[x^i, ChiSquareDistribution[n], x], {i, 
          4}]

Alternately, the moments could be could be obtained from a general result for expected values for positive powers of x.

expect = 
      ExpectedValue[x^i, ChiSquareDistribution[n], x, 
        Assumptions -> {i > 0}]

FullSimplify[expect /. i -> Range[4], 
      Assumptions -> n > 0]

Obtaining the moments from derivatives of the moment generating function is yet another option.

chisqMGF = 
      ExpectedValue[Exp[t*x], ChiSquareDistribution[n], x, 
        Assumptions -> {Abs[t] < 1/2}]

Table[
      Limit[D[chisqMGF, {t, i}] /. t -> k, k -> 0], {i, 
        4}]

Given a random sample of (χ^2) distributed numbers, the observed value of the expectations with the theoretical values.

randchisq = 
        Block[{$RandomState}, SeedRandom[157]; 
          RandomArray[ChiSquareDistribution[20], {10000}]];

Table[ExpectedValue[x^i, randchisq, x], {i, 4}]

evs /. n -> 20

A Polynomial Regression Example

A very common task is to fit a set of data, and then look at some diagnostics to assess the goodness of fit and to check the validity of assumptions about the model.

Here is some data that should fit a quadratic curve in the predictor variable.

noisydata = 
        Block[{$RandomState}, SeedRandom[1234567]; 
          Table[{x, 
              2 + 1.5  x - 3  x^2 + 
                Random[NormalDistribution[0, 2]]}, {x, -2, 
              2,  .1}]];

The Regress function can be used to fit linear models and obtain diagnostics.  Here the best fit curve, standardized residuals, and Cook's distances are returned.

basis = {1, x, x^2};
reg = 
      Regress[noisydata, basis, x, 
        RegressionReport -> {BestFit, StandardizedResiduals, 
            CookD}]

A plot of the data and the best fit curve can be used for a qualitative check of the fit.

<< Graphics`

DisplayTogether[
        ListPlot[noisydata, PlotStyle -> {Red, PointSize[ .02]}], 
        Plot[BestFit /. reg, {x, -2, 2}]];

A Polynomial Regression Example

Other graphics can be used for qualitative checks of fit diagnostics.  For instance, a plot of the standardized residuals can be used to check for trends in the residuals that may indicate a problem with the assumption of constant variance.

ListPlot[resids = (StandardizedResiduals /. reg), 
        PlotStyle -> {PointSize[ .02]}];

A quantile-quantile plot of the standardized residuals against a standard normal distribution can be used for a check of the assumption of normally distributed residuals.

QuantilePlot[resids, 
        Table[Quantile[NormalDistribution[0, 1], y], {y,  .01,  .99,  .01}], 
        SymbolShape -> Point, 
        SymbolStyle -> {PointSize[ .02]}];

A stem plot of Cook's distances can be used to check for points that exert undue influence on the fit, based on some predetermined critical value.

critval = 
      Quantile[FRatioDistribution[Length[basis], 
          Length[noisydata] - Length[basis]],  .5]

MultipleListPlot[CookD /. reg, SymbolShape -> Stem, 
        Epilog -> {Red, Dashing[{ .08,  .03}],  
            Line[{{0, critval}, {Length[noisydata], critval}}]}, 
        PlotRange -> {Automatic, {0, 1}}];

In this case, no Cook's distances are greater than the predetermined critical value.

5.1 Additions and Enhancements

Improved efficiency for data smoothers

Improved efficiency for random number generation from discrete distributions

Extension of CDF for multivariate normal and t distributions to higher dimensions

Addition of clustering functionality

Data Smoothers

Timings are for data of the length given, and length 10 filter windows.  Times are averages of (10^(6 - (Subscript[log, 10])(length))) iterations.

MovingAverage

({{"length", "v 5.0.1", "v 5.1"}, {1000, 0.001094`, 0.00034400000000000007`}, {10000, 0.011410000000000002`, 0.003120000000000001`}, {100000, 0.13899999999999998`, 0.034399999999999986`}, {1000000, 1.2969999999999997`, 0.359`}})

MovingMedian

({{"length", "v 5.0.1", "v 5.1"}, {1000, 0.007906000000000002`, 0.0038280000000000002`}, {10000, 0.1186`, 0.03796999999999999`}, {100000, 1.3672000000000002`, 0.3750000000000002`}, {1000000, 11.311999999999998`, 3.766`}})

Both MovingAverage and MovingMedian are roughly 2 to 3 times faster for the window size chosen.  The speed up is even greater for MovingAverage as the window size increases.

RandomArray for Discrete Distributions

The following timings are for generation of 100,000 random numbers from the specified distributions using RandomArray on a 2.4 GHz Windows machine.

({{"", "v 5.0.1", "v 5.1"}, {BernoulliDistribution[0.5`], 0.141`, 0.046999999999999986`}, {BinomialDistribution[100, 0.5`], 5.922000000000001`, 0.7030000000000001`}, {DiscreteUniformDistribution[100], 0.1559999999999997`, 0.04699999999999993`}, {GeometricDistribution[0.5`], 0.07800000000000029`, 0.07800000000000007`}, {HypergeometricDistribution[50, 40, 100], 140.25`, 0.609`}, {LogSeriesDistribution[0.5`], 17.89099999999999`, 0.25`}, {NegativeBinomialDistribution[100, 0.5`], 13.281000000000006`, 0.8130000000000002`}, {PoissonDistribution[10], 75.344`, 0.5`}})

Timings will vary somewhat with parameter values, but the 5.1 code will be faster than the 5.0 code almost everywhere.

GeometricDistribution shows no speed up because speed enhancements were made to that distribution prior to 5.1.

Multivariate CDFs

Here is a five-dimensional normal example.

Block[{$RandomState}, 
SeedRandom[
          1111]; 
sigma = 
          Table[0, {5}, {5}]; 
Do[
          sigma[[i, i]] = Random[Real, {0, 3}], {i, 5}]; 
        Do[sigma[[i, j]] = 
            Random[Real, {-1, 1}]*sigma[[i, i]]*
              sigma[[j, j]], {i, 5}, {j, 
            i - 1}]; 
sigma = 
          sigma . Transpose[sigma]; 
mu = 
          Table[Random[Real, {-5, 5}], {5}]];

sigma // MatrixForm

mu

Here are 21 cdf values computed for this distribution.  Results tend to be pretty fast for lower PrecisionGoal values.  Sharp features in the integrand near the boundaries of the integral can make computation for higher dimensional integrals at the default PrecisionGoal much slower.

Timing[
      cdfvals = 
        Table[{x, 
            CDF[MultinormalDistribution[mu, sigma], {8, x, 10, 0, 17}, 
              PrecisionGoal -> 3]}, {x, -10, 10, 1}]]

Multivariate CDFs

A ListPlot of the cdf values shows the same shape as a univariate normal cdf as is expected.

ListPlot[cdfvals, PlotJoined -> True];

The last value in cdfvals is effectively the probability that (Subscript[X, 1])<8, (Subscript[X, 3])<10, (Subscript[X, 4])<0, and (Subscript[X, 5])<17, so we can compare it with the cdf of the normal distribution with the second dimension dropped.

cdfvals[[-1, -1]]

CDF[
      MultinormalDistribution[Drop[mu, {2}], 
        Drop[Drop[sigma, None, {2}], {2}]], {8, 10, 0, 17}, 
      PrecisionGoal -> 3]

Multivariate CDFs

For an example of the multivariate t case, we can convert the covariance matrix into a correlation matrix.

MatrixForm[
      corrmat = 
        Transpose[sigma/Sqrt[Table[sigma[[i, i]], {i, 5}]]]/
          Sqrt[Table[sigma[[i, i]], {i, 5}]]]

Timing[
      CDF[MultivariateTDistribution[corrmat, 100], {8, 2, 10, 0, 17}, 
        PrecisionGoal -> 3]]

For equivalent dimensions, the t integral is more complicated than the normal integral and will be slower in general.  Since the degrees of freedom are fairly large in this example, the result could be approximated using a multivariate normal cdf.

Timing[
      CDF[MultinormalDistribution[{0, 0, 0, 0, 0}, corrmat], {8, 2, 10, 0, 
          17}, PrecisionGoal -> 3]]

Clustering

Clustering functionality is included in the Statistics`ClusterAnalysis` standard add-on package in Mathematica 5.1.  This includes clustering of numeric, Boolean, and string data.

Here are ten sets of bivariate normal data points.

origclusters = 
        Block[{$RandomState, muvals, sigmavals}, 
          SeedRandom[1234546789]; 
muvals = 
            Table[Random[
                Real, {-10, 
                  10}], {10}, {2}]; 
sigmavals = 
            Table[Block[{s1 = Random[], s2 = Random[], 
                  s12}, 
s12 = 
                  Random[]*Sqrt[s1*s2]; 
{{s1, s12}, {s12, 
                    s2}}], {10}]; 
Join[
            MapThread[
              RandomArray[MultinormalDistribution[#1, #2], 
                  Random[Integer, {20, 75}]] &, {muvals, 
                sigmavals}]]];

This function will differentiate lists of points by color.

ClusterPlot[cl_, opts___] := 
        Module[{colors, len = Length[cl]}, 
colors = 
            Map[Hue, Range[1, len]/(len + 1)]; 
Show[
            Graphics[{AbsolutePointSize[5], 
                Transpose[{colors, Map[Thread[Point /@ #] &, cl]}]}, opts,  
              Frame -> True]]
];

Here are the original clusters.

ClusterPlot[origclusters];

Clustering

This will combine the points in the original clusters into one list of points and reorder them.

Off[General::shdw]; Needs["DiscreteMath`Combinatorica`"]

data = 
        Block[{$RandomState}, RandomSeed[1]; 
          RandomPermutation[Flatten[origclusters, 1]]];

Here we can see the clusters found by FindClusters.

ClusterPlot[FindClusters[data]];

FindClusters tries to find an optimal number of clusters by default.  Alternately, the desired number of clusters can be specified.

Here are the clusters found if we require exactly five clusters.

ClusterPlot[FindClusters[data, 5]];

Clustering

The Agglomerate function uses agglomerative hierarchical clustering to obtain a full cluster hierarchy.  DendrogramPlot gives a graphical representation of that hierarchy.

Here is the first paragraph from section 0.2.1 of The Mathematica Book split into individual strings.

mbooktext = 
      StringSplit["Mathematica is the world's only fully integrated \
environment for technical computing. First released in 1988, it has had a \
profound effect on the way computers are used in many technical and other \
fields."]

Here is a graphical representation of the hierarchical clustering of those strings.

DendrogramPlot[mbooktext, LeafLabels -> (# &), 
        Orientation -> Left, ImageSize -> 600];

Summary

We have seen a very brief introduction to statistical functionality built into Mathematica.  We have also seen new statistical features added in 5.1, as well as improvements to functions that existed prior to 5.1.

For those interested in a more complete coverage of statistics in Mathematica, the Wolfram Education Group currently offers a mini-course on Applied Statistics with Mathematica specific to regression, which is currently offered online.  An Applied Statistics with Mathematica mini-course on descriptive and mathematical statistics is also in development.  There are also a number of books about statistics with Mathematica available.