mathStatica and Assumptions Technology

Colin Rose and Murray D. Smith

Abstract
The mathStatica pacakage provides a completely general toolset for doing mathematical statistics with Mathematica. Common statistical operators are defined for taking expectations, finding probabilities, deriving transformations of random variables and so on... all from first principles. By taking full advantage of the latest Version 4 Assumptions technology, these functions automatically yield exceptionally clean and neat symbolic output.

Dr Colin Rose  	  Theoretical Research Institute,  Sydney	colin@tri.org.au
Dr Murray D. Smith	  University of Sydney  			Murray.Smith@econ.usyd.edu.au

Introduction

.	A New Approach

The use of computer software in statistics is far from knew. Indeed, literally hundreds of statistical computer programs exist. Yet, underlying existing programs is almost always a numerical / graphical view of the world. Mathematica can easily handle the numerical and graphical sides, but it offers in addition an extremely powerful and flexible symbolic computer algebra system. The mathStatica software package builds upon that symbolic engine to create a sophisticated toolset specially designed for doing mathematical statistics. This forms the basis of a Springer-Verlag text entitled Mathematical Statistics with Mathematica (600 pages) that ships with CD-ROM, custom palettes, online help, and live interactive chapters. This provides a simple and unified approach for doing mathematical statistics with Mathematica, suited to those in statistics, econometrics, engineering, physics, psychometrics, economics, finance, biometrics - indeed across the full ambit of the physical/social sciences and across both the pure and the applied domains. It is expected to be available by January 2000.

.2	Features

Features include: * a complete suite of functions for manipulating probability density functions, * symbolic maximum likelihood estimation, * numerical maximum likelihood estimation, * automated symbolic Pearson curve fitting, * Johnson fitting, * Gram-Charlier expansions, * nonparametric kernel density estimation, * moment conversion formulae (convert automatically between cumulants, raw moments, central moments, and factorial moments: univariate and multivariate), * random number generation for any discrete distribution, * fully automated transformations (functions of random variables), * asymptotics, * decision theory (order statistics, mean square error, ...), * unbiased estimation (h-statistics, k-statistics, polykays). The mathStatica package can replace many reference works and extends the analysis to problems of arbitrary high order. This paper concentrates on just 4 functions, and illustrates how Assumptions technology is used to automatically generate crisp textbook quality answers.

.3	Design Philosophy

The mathStatica package has been designed with two central goals:  it sets out to be general, and it strives to be 'delightfully simple'.

2	Core Functions

2.1	Getting started

The mathStatica package adds over a hundred new functions to Mathematica. But 95% of the time, we can get by with just 4 of them:

({{"\!\(\*StyleBox[\"\\\"abbreviation\\\"\", \"SO10\"]\) ", "\!\(\*StyleBox[\"\\\"description\\\"\", \"SO10\"]\) "}, {PlotDensity[ f ]       , Plotting (automated)}, {Expect[x,f ] , Expectation operator        E[X] }, {Prob[x,f], \(Probability\ \ \ \ \ \ \ 	\ \ \ \ \ \ \(P(X ≤ x)\)\)}, {Transform[eqn,f ]      , Transformations}})

This ability to handle plotting, expectations, probability, and transformations, with just 4 functions, makes the mathStatica system very easy to use, even for those not familiar with Mathematica.

To illustrate, let us suppose the continuous random variable X has probability density function (pdf) f(x)  =   (1/(π (Sqrt[ 1 - x]) (Sqrt[x]))), where x  ∈  (0, 1). In Mathematica, we enter this as:

f  = 
      (1/(π (Sqrt[1 - x]) (Sqrt[x])));            domain[
        f]  =  {x,  0,  1};

This is known as the Arc-Sine distribution. Here is a plot of f(x):

PlotDensity[f];

Here is the cumulative distribution function (cdf), P(X ≤ x), which also provides the clue to the naming of this distribution:

Prob[x,  f]

The mean, E[X], is:

Expect[x,  f]

while the variance of X is:

Var[x,  f]

The (r^th) moment of X is E[(X^r)]:

Expect[(x^r),  f]

The moment generating function (mgf) of X is E [(e^(t X))] :

Expect[(E^(t x)),  f]

g   =  
      Transform[y  == (Sqrt[x]),  f]
domain[g]  =  TransformExtremum[y  == (Sqrt[x]),  f]

So, we have started out with a quite arbitrary pdf f(x), transformed it to a new one g(y), and since both density g and its domain have been inputted into Mathematica, we can also apply the mathStatica tool set to density g.

2.2	Introducing Assumptions Technology

The mathStatica package has been designed to seamlessly support parameters. It does so by taking full advantage of the new Assumptions technology introduced in Version 4 of Mathematica. To illustrate, let us consider the familiar Normal distribution with mean μ and variance (σ^2). That is, let X~N(μ,  (σ^2)), where μ  ∈  ℝ and σ > 0. We enter the pdf f(x) in the standard way, but this time we have some extra information about the parameters μ and σ which can be added to the end of the domain[f] statement:

f    =    (1/(σ (Sqrt[ 2 π]))) 
          Exp[-((((x  -  μ))^2)/(2 (σ^2)))];      
domain[
          f]  =  {x,  -∞,  ∞}  &&  {μ ∈ Reals,  σ > 0};

From now on, the assumptions about μ and σ will be 'attached' to density f, so that whenever we operate on density f with a mathStatica package function, these assumptions will be automatically applied in the background. With this new technology, the mathStatica package can usually produce remarkably crisp textbook-style answers, even when working with very complicated distributions.

The mathStatica package function, PlotDensity, makes it extremely easy to examine the effect of changing parameter values. For instance, here is the pdf plotted at three different values of σ:

PlotDensity[
        f  /.  {μ -> 
              0,  σ ->  {1,  2,  3}}];

It is well known that E[X]  =  μ and Var(X) =  (σ^2), as we can easily verify:

Expect[x,  f]

Var[x,  f]

Because the mathStatica package is completely general in its design, we can just as easily solve problems that are both less well-known and more 'difficult', such as finding Var((X^2)) :

Var[(x^2),  f]

or find the entropy, given by E[-log(f(x))]:

Expect[-Log[f],  f]

Assumptions technology in mathStatica:

1. mathStatica implements assumptions technology in a distribution-specific manner. That is, the assumptions are attached to the density f(x;  θ), in contrast to the parameter  θ itself.

2. This means we can work with two or more pdf's, say:

	f(x;   θ) 	assuming parameter θ > 0,      and simultaneosuly 
	
	h(x;   θ) 	assuming parameter θ < 0 

and the package will automatically assume that θ is positive when working with f(x), but negative when working with h(x).

3. Finally, because the assumptions are not global, it is important to realise that the assumptions will only be automatically invoked when using the suite of mathStatica package functions. By contrast, Mathematica's inbuilt functions such as the derivative function D[f, x], will not automatically assume that θ  > 0 when presenting the solution. ■

2.3	Discrete random variables

The mathStatica toolset automatically handles discrete random variables in the same way. The only difference is that, when we define the density, we add a flag to tell Mathematica that the random variable is {Discrete}. For instance:

Let the discrete random variable X have probability mass function (pmf) f(x) =  (({{r + x - 1}, {x}})) (p^r) (((1 - p))^x), where x = 0,  1,  2,   ... . Here parameter p  is the probability of success, while parameter r is a positive integer. In Mathematica, we enter this as:

f   =   Binomial[r + x - 1, 
            x] (p^r) (((1 - p))^x);
domain[
          f]  =  {x, 
            0, ∞}   &&  {Discrete}   && 
                                                  {0 < p < 1,   r  > 0,  
            r ∈ Integers} ;

This is known as the Pascal distribution. Here is a plot of f(x):

PlotDensity[
        f /.  {p ->  (1/2),  r ->  10 }];

Here, for example, is the cdf, equal to P(X ≤ x),

Prob[x,  f]

The mean E[X] and variance of X are given by:

Expect[x,  f]

Var[x,  f]

The probability generating function (pgf) is E [(t^X)] :

Expect[(t^x),  f]

2.4	Multivariate pdf's

The mathStatica toolset extends naturally to a multivariate setting. To illustrate, let us suppose that X and Y have joint pdf f(x,  y) with support x > 0, y > 0:

f    =   (E^(-2 (x + y))) ((E^(x + y)) + (-2 + (E^x)) (-2 + (E^y)) α);
domain[
          f]  =  {{x,  0,  ∞},  {y,  
              0,  ∞}}        &&     {-1  <  α < 1};

where parameter α is such that -1 < α < 1. This is known as Morgenstern's (1956) Bivariate Exponential pdf. Here is a plot of f(x,  y):

Plot3D[f /.  α ->  - .8, {x,  0,  2}, {y,  0,  1}, 
AxesLabel -> {"\!\(\*
StyleBox[\"x\",\nFontSlant->\"Italic\"]\)", "\!\(\*
StyleBox[\"y\",\nFontSlant->\"Italic\"]\)", "\!\(\*
StyleBox[\"f\",\nFontSlant->\"Italic\"]\)"}];

Here is the cdf, namely P(X ≤ x,  Y ≤ y),

Prob[{x,  y},  f]

Here is Cov(X,  Y), the covariance between X and Y:

Cov[{x,  y},  f]

More generally, here is the variance-covariance matrix:

Varcov[f]  //  MatrixForm

Here is the marginal density of X:

Marginal[x,  f]

Here is the conditional pdf of Y, given X = x:

Conditional[x,  f]

Here is the bivariate mgf \(TraditionalForm\`E[e\^\(\(t\_1\) X\  + \ \(t\_2\) Y\)]\):

mgf  =  
      Expect[ (E^((Subscript[t, 1]) x  +  (Subscript[t, 2]) y)),  
        f]

Differentiating the mgf is one way to derive moments. Here is the product moment E[(X^2) (Y^2)]:

\(D[mgf, \ {t\_1, 2}, \ {t\_2, \ 2}] /. \ t\__ -> \ 0\  // \ 
      Simplify\)

which we could otherwise have found directly with:

Expect[(x^2) (y^2),  f]

eqn  =  {u ==  (y/(1 + x)),  
          v  == (1/(1 + x))};

Using Transform, we can find the joint density of random variables U and V, denoted g(u,  v):

g  =  Transform[eqn,  f]

while the extremities of the new random variables are:

TransformExtremum[eqn,  f]

2.5	Bipartite pdf's  (kinks)

Some density functions take a bipartite form. To illustrate, let us suppose X is a continuous random variable x  ∈  (0, 1) with pdf:

where 0 < α < 1. We enter this as:

f     =     If[
          x < α , 2 (( α - x)/α), 
          2 (( x - α)/(1 - α))] ;   
domain[
          f]   =    {x, 0, 
            1}   &&   {0 < α < 1};

This is known as the Inverse Triangular distribution, as is clear from a plot of f(x):

PlotDensity[
        f /.  α -> {(1/4), (1/2),  (3/4)}];

Here is the cdf, P(X ≤ x):

Prob[x,  f]

Note that the solution depends on whether x < α or x ≥ α. Here is a plot of the cdf at the same three values of α :

cdf =  Prob[x,  f];
Plot[{cdf /. α -> (1/4),  
            cdf /. α -> (2/4),  cdf /. α -> (3/4)}, {x,  0,  1},AxesLabel -> {"\!\(\*
StyleBox[\"x\",\nFontSlant->\"Italic\"]\)", "\!\(\*
StyleBox[\"F\",\nFontSlant->\"Italic\"]\)"}, PlotRange ->  All, PlotStyle  ->  {AbsoluteThickness[2],  {},  
              AbsoluteDashing[{7, 5}]}];

The mathStatica toolset still operates in the standard way. For instance, the mean E[X] and variance are given by:

Expect[x,  f]

Var[x,  f]

and here is the entropy, given by E[-log(f(x))]:

Expect[-Log[f],  f]

3	Acknowledgements

The mathStatica package was partly developed while Rose was visiting Wolfram Research, and while Smith was visiting the Sonderforschungsbereich 386, Institüt für Statistik, Ludwig-Maximilians-Universität, München. Rose gratefully acknowledges a Visiting Scholar Grant from Wolfram Research, while Smith acknowledges gratefully the support of the Alexander von Humboldt-Stiftung, and the Australian Research Council (through VISLAB and MEMLab). ■