NIACSYNIA as a Tool for Statisticians

J. Michael Steele, Princeton University

ABSTRACI this report is an artifact of using the typesetting features

Three concrete examples are used to show how sym of MACSYMA ( check out the flag setting typeset:true).

bolic: computation can be of value in statistical research. (cl) difF(sin(x)/xx,3~,

Enough of the syntax of MACSYMA is discussed to pro~

vide an introduction to actual usage. The questions con

sidered concern the relationship of conditional cumulants 3 sin(x 6 sin(x cos(x ) + 6 cos(x (d 1)

Of homogeneous statistics, a method for proving limit X X X

theorems by studying zeros of polynomials, and the

Inendre transformation of some probability densities. (c2) taylor(cos(x)^ 3,xAQ

One question which is constantly considered is * When

will MACSYMA seriously help 3x 2 + 7x 4 _ 61xl + (d2)

1. Introduction 2 8

MACSYMA and other symbolic computational tools

are now widely available. Still, the statistical commun (c3) solve(2*x^ 2+x+2Y

ity has not felt the full impact of these new tools and

many people may not yet realize that the conduct of '115i + 1 J i~i 1

research in probability and statistics will be forever [X ,X 1 WM

changed by their presence. The purpose of this report is 4 4

to illustrate how even naive use of MACSYMA can be of

substantial benefit in research activity which is far from (cS) quit(Y,

being a priori computational.

The plan of the report is to first give a gentile intro

duction to the structure and syntax of MACSYMA, give This microsession is enough to show some of the

a few hints for using MACSYMA comfortably In a power in symbolic computation. One does not need

[INIX environment, and then give three examples. The MACSYMA to calculate the third derivative of sin(x)/x,

are not concocted but are just selected from the or the Taylor series Of COS(X )21 to four terms, but it is in

rny own experience in using MACSYMA as one of the fact easier and more reliable than doing it by hand.

tools of daytoday research. Finally, brief comments are Definite and indefinite integrals are also available with a

given con(erning the ways one can customize a comparably natural synt,m The use of solve illustrates a

\4;\('SYMA environment to facilitate the types of com slight variation on natural syntax but

putdtJons common in probability and statistics. solve(2*x 2+x+3)=0 works just as well.

Before going to a real example, one should learn

2. Getting Started some quick syntax. First "*,* terminates a command (or c

The way to get started using MACSYMA is to fire ) line. Also, [a,bA is a list containing a, b, and C. There

are assignment operators *." which defines functions

up MACSYMA, open An Introduction to MACSYMA

by Joel Moses and the MACSYMA Group of Synibloics~ and for more typical assignment.

Ind1984), and start playing. This is fun and effective. In 3. Some Practical Hints

!stark contrast, you can go to MACSYMA with a tough There are some frustrations which most newcomers

problem which you have not succeeded in doing by

hand, try to use MACSYMA for the first time and end to MACSYMA seem to meet and which are easily over

up meeting only frustration and likely failure. come once 'you know the trick*. First, the beginner

should avoid the N4ACSYMA editor. The best way to use

In teaching MACSYMA to my friends, the hardest MACS~ Is to work interactively only while explor

thing to get across always seems to be that MACS~ is

!rig and playing. Once a serious problem is to be studied

just a tool. A proper attitude Is created by thinking of the MACSYMA commands should be put in a file, say

MACSYMA as being like integration by parts. Surely

compfile, and executed in MACS~ with the com

integration by parts is a magnificent tool. It shines in mand batchCcomphie)

vital areas like the calculus of variations. It gives birth to The benefits of this process are (1) If an error is

the formulas of Gauss, Green, and Stokes. It is a weapon

made you can stop MACS~ with ^ Z and then just

4 first enoill in many parts of asymptoties. But. for all edit compfile width your favorite editor, (2) The code

of thw virlues, integration by parts will not solve dil you create becomes an asset which can be used in

prnblexiis, and neither w ill MAC.SYMA. The first task of

difrerent contexts after the traditional pattern of

the beginner is to shitn the inevitable disappointment of "modify rather than writew, and (3) The code can be

naive P ~*~(tations and to 4rticulate those areas of per conveniently quoted In papers, shared with others, etc.

sonal mk.erest where MACS1i'MA has a. realistic chance of

Please do not be mislead. There are many other

tools in MACSYMA which allow one to communicate

Fhewope of MA(SYMA and its interactive nature with files and with structures within MACSYMA. but.

are casv to see by looking over a brief but complete

for a while, these seem best left to vertuosi.

UNIX script of a trial session. In the first computations

we review, the syntax should be obvious if one bears in Two more practical hints ? Even if you do not wish

mind that the clines are entered by the user and the d to do ANY symbolic computation MACSYNIA can be a

lines are MACSYNIA's responses. One small point is that ('y~nd for Its FORTRAN and eqn facilities. in a nut

^ is used to denote exponentiation. The symbols W11 Itell 1.()R'rRAN(exi)ression) produces the FOR1RAN

W2), etc. will appear on the left side on the terminal just (ofle it) compute the specified expression. 1 fin cl this a

like the clabels. Their appearance on the right side in v(r~ pleasant way to avoid coding mistake., In the same

way the typeset flag of MACSYNIA permits, you to get 'he

t e.qn code of the resulting dline. This is much more pleasant and precise than writing your own equations with eqn. The equations of this report were all done in this way.

4. Simplicity 'rhough Diligence

One way in which MACSYNIA can be of aid in research is simply by removing the difficulties of looking at the problem from many sides. Often one has an idea which can go unexplored because the calculations look too cumbersome. This inhibition emerges for two reasons. First, it takes conviction to go trough a tedious calculation which may not lead to a useful result and second, the notion of usefulness is often linked to simplicity.

This block can be busted by MACSYMA in a number of cases. One from my experience concerns conditional cumulants of homogeneous statistics. The reason on can expect MACSYMA to be efFective here is that the tasks required are just tedious, not ingenious. At difFerentiation, substitution, and rational simplification by expansion MACSYMA is as surefooted as a mountain goat.

5. Cumulant Correction Formulas

We will assume that f satisfies the homogeneity condition (1.1) and set

Z', f(X1X2,,M (2.1)

where the Xi are independent and uniform on [0,11d. Further, we denote by Z, a random variable with the same definition as Z,, except that one of the Xi's is constrained to be on "forward" part of the boundary of [0, 1 ]d; that is, one of the Xi's is constrained to be in the Set B _ [0, Ild [0, 1)d.

If Y is any random variable the cumulant generating function of Y is defined by

7_Ki(Y) L log EW') (12)

i =1 i !

and the numbers Kl(Y) are the cumulants (or semiinvariants) of Y. Naturally, Kl(Y) is the mean Of Y, KZ is the variance, and these are the most important of the cumulants. The third and forth cumulants ( the skewness and the kurtosis) are of interest principally because after the variance all of the cumulants of a normal distribution are zero. The cumulants of Z, do not have to tend to zero for it to be asymptotically normal, but if the third and forth cumulants tend to zero it is a strong in(ii(.dlion of asymptotic normality.

One of our goals in pursuing the simulation method stimlied here (and in particular with bothering about higlier cumulants) is to obtain a compulationally ,.\1)~;dieni rnethod of assessing the possibility of asymptot ic noi inality in the case of computationally important random variables like the length of the minimal spanning trm In Rainey (1982) the asymptotic normality of the MST was proved under the hypothesis of an assumplion which seems very difficult to verify (but which is still quite plausible). The asymptotic normality of the length of the MST remains an important open problem. In fact, there is no subadditive Euclidean functional for which normality has been roved yet it is quite likely to hold for the whole class (` c.f. Beardwood, Halton, and flammersley (1959) and Steele (1981)).

To begin the conditioning argument, we let

S = minis: Xi E [O,Sld) thus, S is the size of the side of the smallest cube contain

i

i

i

1

2

!rig the data and the origin. We note that with probability one there is exactly one of the Xi on the forward part of the boundary of [O,Sld. A key observation is that we then have that the distribution of Z,, equals the distribution of SPZ,, where 5 and Z,, are generated independently.

Now we can just calculate.

Z,

exp K, (Z" ) 1, T ~ E(EW ~S)) (2.3)

Me"'Z"IS)) EexpFKi(Z,,)S'PL~

i=l

Next we expand both sides of the above using

g 1 2

exp F a, t + a, + (a2 + a2)!,'

i =1 7! Z

+ (a 3 + 3a2a l + a3)iii +

1 3!

(a 4 + 6a2a 2 + 4aa, + 3a 2 tAl

1 1 2 + a,) ~! +

and then complete taking the expectation by using the elementary fact that E (SO') dn (ot + dn )1. On simplifying the algebra, one finds the following tidy expressions for the conditional cumulants (To save space we will write Ki for Ki (Z, ) :

KI(Z',,) = K10 + L) (2.4a)

nd

K2(Z', K2(1 + 2P) K 2 p (2.4b)

nd n 2 d

K3(7.'n) = K3(1 +

and finally we have

3p 6K 1K2p2 2K3 3

+ 1 p 2.4c)

nd nTdT illTI71

1+ 4p

dn

2 3 4 4

+ 24K1K2P 6KIP

n d

2) P 2

12(K1K3 + KZ F 2.4d) n rd~

The method we have used is a natural one, but there are two surprises. The first of these is that the homogeneity hypothesis is powerful enough to give such simple results. The second is that it turns out to be incredibly neater to give an expression for the conditional variances in terms of the unconditional ones rather than viceversa.

The idea of using a generating function approach to obtain conditional cumulant formulas is inherent in the definition of cumulants. It has been used often, and in particular the method is applied in Brillinger(1969) to obtain a general formula relating cumulant to conditional cumulants. Interestingly, that an attempt to obtain the identities (2.4ad) directly from Brillinger's general formula is almost certainly doomed to failure because of the very much greater algebraic complexity of the expressions for the unconditional cumulants in terms of the conditional ones.

rhere is a second method for obtaining results relating conditional and inconditional cumulants which uses Mobius function of the lattice of partitions. Speed(I983) introduced this method and showed its attractiveness in a number of examples including an intriguing derivation

of Brillinger's formula. The Moblus function methoa For example, the following famous transform relation is

shows great promise; but, at the present level 1 of complex the key to Holders inequality,

ity, series inversion used here still seems to be the tool of

choice. XY < _XP + Y~

Some comments on the formulas (2.4ad) are in p q

order before going on to the simulation method. In the provided 1+ 1. In the same way that one can use

first place, even the modest equation (2.4b) is able to pro p q

vide some interesting lower bounds on the variance of Holder's inequality to determine the the dual space of

the pth power integrable functions, one can start off

certain subadditive Euclidean functionals. We will illus with any Legendre transform pair.

trate with the traveling salesman problem (TSP) in

dimension 2 , but comparable results will obviously hold A statistical problem of interest is the determination

for the minimal spanning tree, minimal triangulation, of all of those functions h (g,o) which have a finite

and other functionals and other dimensions. integral when multiplied by the normal density

The only fact we require about the TSP is that if (b(x ;g, a) and integrated with respect to g and a. To

(X lly 2, '. ' X, ) is the length of the shortest tour stud this MACSYMA was used to calculate the (bivari

through the points Xi 1 Q ~i then K1 is asymptotic to ate ~ Legendre transformation of (t with respect to

1 A and a. The result takes a page to express. Still, this

~n _f as n tends to infinity (Beardwood, Halton, Ham worm of an expression is probably the best available tool

for determining which priors on p and a will guarantee

mersley (1959)). Since this f is homogeneous of order

1 a posterior distribution which is absolutely continuous

p=l, equation (2.4b) implies that K2 is fX ). This modest with respect to Lebesgue measure.

n

bound is of interest because it is known that K2 IS For a discussion of the classical Legendre transfor

[x)unded for all n ( Steele (1981) ), but no better lower mation ( as it applies to Hamiltonian systems ) as well as

bound than the one just give is known. Naturally, it is one version of a MACSYMA Legendre transformation

expected for the TSP in d=2 that K2 converges to a posi one should consult Rand (1984).

tive constant as n goes to infinity. Mixed Computations

6. Simplicity Though Diligence The symbolic computational capacities Of

MACSYMA are augmented with numerical and graphi

One way in which MACSYMA can be of aid in cal capabilities. This addition to the workbench permits

research is simply by removing the difficulties of looking

the joint investigation of "formulas" and the exploration

at the problem from many sides. Often one has an idea of their quantitative features.

which can go unexplored because the calculations look

too cumbersome. This inhibition emerges for two reasons. One way this plays out to the investigators advan

First, it takes conviction to go trough a tedious calcula tage is that it is possible to gain insight from much larger

tion which may not lead to a useful result, and second, examples than would normally be possible. A useful

example of this from my own experience come for the

the notion of usefulness is often linked to simplicity. theory of random trees.

This block can be busted by MACSYMA in a The detailed example is given in Steele(1985), but

number of cases. One from my experience concerns con the role of MACSYMA is easy to review without those

ditional cumulants of homogeneous statistics. The reason details. The generating function for the number of leaves

on can expect MACSYMA to be efective here is that the of a random tree was derived in Renyi(I959) under the

tasks required are just tedious, not ingenious. At model of randomness which corresponds to the uniform

differentiation, substitution, and rational simplification

distribution on the Pruffer codes ( These codes provide a

by expansion MACSYMA is as surefooted as a mountain natural onetoone correspondence between labeled

goat.

unrooted trees and n2 tuples of the integers 1 <, k ~n

One of the properties of mathematics ( which also For the purpose of modeling, the uniform model on

persists in statistics) is a strong preference for simple

trees is quite limited , eg. one can not escape the conse

~ers. Part of the process ot getting to understand the

n

power of symbolic computation rest in comming to terms quence of there being approximately leaves on the

e

with complex results. What makes this novel is that one great majority of the trees. flow much nicer for model

typically does not have to deal with complex results ing purposes if we had a family of tree valued random

because once a calculation gets past a certain complexity

variables whose expected number of trees could be

threshold it is given up for lost. In a computational chosen by the modeler ? There is such a family and a

environment one can get extremely complex results and method for generating choices from that family given in

one faces the new challenge of using such results intelli

gently. Steele(1985).

One virtue of Renyi's distribution on random trees

is that Reny! was able to prove a central limit theorem

the number of leaves of a tree form his family ( as n

7. The Legendre Experience tends to infinity). To establish the corresponding result

in the exponential family of random trees it appeared to

This feature of the changing role of complexity is be almost necessary to show that the generating func

well illustrated by the enjoyable practice of computing tions given by Renyi had all real zeros. Naturally that

Inendre transformations. For any function f(x) the opinion takes considerable motivation. Using MACSYNIA

Lnendre transformation may be defined by it was possible to examine polynomials up to n=10 were

g (y max (xy f (x D it was discovered that for large n there were nontrivial

complex roots. The plan of analysis was saved when it

The importance of this transformation in analysis was found that by adding a stochasticly negligible term

comes In part from the relaxed inequality

to Renyi's polynomial the roots could be forced back to

ry < f (X) + g (y) the real line. This eventually lead to a proof of the cen

3

tral limit theorem for the number of leaves of a tree Gladd, N. T. (1984) Symbolic Manipulation Techniques

from the exponential family. This is particularly satisfy for Plasma Kinetic Theory Derivations Journal of Com

ing because it plays back into the original motivation by Putational Physics, 56, 175202.

making the mean number of leaves an essential parame

ter. Hochbaum, D. and Steele J. M. (19821 Steinhaus's

Q. Concluding Remarks geometric location problem for random samples In the

1 This introduction to MACSYMA has been casual plane, Adv. App4 Prob, 14, 5667.

and discursive. The goal has been to motivate others to Karp, R. M. (1977). Probabalistic analysis of partitioning

begin productive interaction with the tools which are algorithms for the traveling salesman problem in the

now widely available. The hints which have been given plane. Math. Oper. Res. 2 202224.

are those that seem to be responsible and welcoming.

There is a great deal of benefit to the statistical commun MACS~ Group of Symbolics, IncPrincipal writer

ity which can be gained from symbolic computation and Joel Moses(I984), An Introduction to MACSYMA Sym

it is inevitable that current elforts will be made to look bolics. Inc. Cambridge Mass.

primative by those which will follow very shortly. It Is

especially intersting to imagine how our collective con

ception of what constitutes a useful answer will migrate Mathlab Group (1983) MACSYMA Reference Manual,

in the next few years. version Ten, Vol. I & 11, Laboratory for Computer Science

ACKNOWLEDGEMENT MIT, Cambridge~ Man

This research was supported in part by NSF grant

DMS8414069. Results reproduced here were obtained Ramey. D. (1982) A Nonparametric Test of Bimodality

in part by using MACSYMA, a large symbolic manipula with Application to Cluster Analysis. Department of

tion program developed at the MIT laboratory for com StatisticsYale University, New Haven Conn.

puter Science and supported from 1975 to 1983 by Rand, R. FL (1984) Computer Algebra in Applied

NASA under grant NSG 1323. by ONR under grant

N0001477C0641. by DOE under grant ET78C02 Mathematics. an Introduction to MACSYMA Pitman

4687,by the USAF under F4%2079C020. and since Advanced Publishing Program, Boston.

1982 by Symbolics, Inc. of Cambridge, Mass. MACSYMA

is a trademark of Symbolics, Inc. Key Words. Reny~ A. and Sulanke, RX19641 Uber die konvexe Hulle

MACSYMA, symbolic computation, cumulants, Legendre von n zufallig gewahIten Punkten. IL Z Wahrscheinli

transformation, central limit theorem. AMS Subject chtkeitsteorie, 3 138147

Classifi cations (19801 Major 68C20; Minor 62EX

Author's Address: Program in Engineering Statistics . Renyi, A~1956). Some Remarks on the Theory of Trees.

Fngineering Quadrangle, Princeton University, Princeton Pub. Math. Inst. Hung. Acad. Sci. 4, 7385.

N. J. 08544 Speed, T. (1983). Cumulants and Partition Lattices. Aust.

J. Statist. 25, 378388

REFERENCES Steele, J. NI. (1980). Subaditive Euclidean Functionals

and nonlinear growth in geometric probability. Ann.

Beardwood, J. Halton, H. J. and Harnmersley, 1M. (1959). Probability. 2 7585.

The shortest path through many ponts. Proc. Camb.

Phil. Soc. 55, 299327. Steele, J. M. (1985). Exponential Family of Trees. Techni

cal Report, Program in Engineering Statistics and

Brillinger, D.R. (1969). The calculation of cumulants via Management Science, Princeton University.

conditioning. Ann. Inst. Statist. Math., 21, 375390.

4

i

1 1