XLISP-STAT is a statistical environment built on top of the XLISP programming language. This document is intended to be a tutorial introduction to the basics of XLISP-STAT. It is written primarily for the Apple Macintosh version, but most of the material applies to other versions as well; some points where other versions differ are outlined in an appendix. The first three sections contain the information you will need to do elementary statistical calculations and plotting. The fourth section introduces some additional methods for generating and modifying data. The fifth section describes some additional features of the Macintosh user interface that may be helpful. The remaining sections deal with more advanced topics, such as interactive plots, regression models, and writing your own functions. All sections are organized around examples, and most contain some suggested exercises for the reader.
\piThis document is not intended to be a complete manual. However, documentation for many of the commands available is given in the appendix. Brief help messages for these and other commands are available through the interactive help facility described in Section 5.1 below.
XLISP itself is a high level programming language developed by David Betz and made available for unrestricted, non-commercial use. It is a dialect of Lisp, most closely related to the Common Lisp dialect. XLISP also contains some extensions to Lisp to support object-oriented programming. These facilities have been modified in XLISP-STAT to implement the screen menus, plots and regression models. Several excellent books on Common Lisp are available. One example is Winston and Horn [22]. A book on XLISP itself has recently been published. Unfortunately it is based on XLISP 1.7, which differs significantly from XLISP 2.0, the basis of XLISP-STAT 2.0.
Documentation on XLISP 2.0 is included in the source code \pndistribution; a copy is available from me. Betz has also written a tutorial introduction to XLISP for BYTE [5]. However, this deals with version 1.2, which again differs significantly from version 2.0.
\piXLISP-STAT was originally developed for the Apple Macintosh. It is now also available for UNIX systems using the X11 window system, for Sun workstations under the SunVIew window system, and, with only rudimentary graphics, for generic 4.[23]BSD UNIX systems. The Macintosh version of XLISP-STAT was developed and compiled using the Lightspeed C compiler from Think Technologies, Inc. The Macintosh user interface is based on Paul DuBois' TransSkel and TransEdit libraries. Some of the linear algebra and probability functions are based on code given in Press, Flannery, Teukolsky and Vetterling [14]. Regression computations are carried out using the sweep algorithm as described in Weisberg [21].
This tutorial has borrowed several ideas from Gary Oehlert's \IMacAnova User's Guide\i [13]. Many of the on-line help entries have been adopted directly or with minor modifications from the \IKyoto Common Lisp System\i. Most of the examples used in this tutorial have been taken from Devore and Peck [11]. Many of the functions added to XLISP-STAT were motivated by similar functions in the \IS\i statistical environment [2,3].
The present version of XLISP-STAT 2.0 seems to run fairly comfortably on a Mac II or Mac Plus with 2MB of memory, but is a bit cramped with only 1MB. It will not run in less than 1Mb of memory. The program will occasionally bomb with an ID=28 if it gets into a recursion that is too deep for the Macintosh stack to handle. On a 1MB Mac it may also bomb with an ID=15 if too much memory has been used for the segment loader to be able to bring in a required code segment.\pn
\piDevelopment of XLISP-STAT was supported in part by grants of an Apple Macintosh Plus computer and hard disk and a Macintosh II computer from the MinneMac Project at the University of Minnesota, by a single \pnquarter leave granted to the author by the University of Minnesota, by grant DMS-8705646 from the National Science Foundation, and by a research contract with Bell Communications Research..
\B\fs<14>Using this Tutorial\b\fs<12>
The best way to learn about a new computer program is usually to use it. You will get most out of this tutorial if you read it at your computer and work through the examples yourself. To make this easier the named data sets used in this tutorial have been stored on the file \B\ff<(CG)Courier>tutorial.lsp\b\ff<(CG)Times> in the \BData\b folder of the Macintosh distribution disk. To read in this file select the \BLoad\b item on the \BFile\b menu. This will bring up an \BOpen File\b dialog window. Use this dialog to open the \BData\b folder on the distribution disk. Now select the file \B\ff<(CG)Courier>tutorial.lsp\b\ff<(CG)Times> and press the \BOpen\b button. The file will be loaded and some variables will be defined for you.\SP1\Sp
\B\fs<14>Why XLISP-STAT Exists\b\fs<12>
There are three primary reasons behind my decision to produce the XLISP-STAT environment. The first is to provide a vehicle for experimenting with dynamic graphics and for using dynamic graphics in instruction. Second, I wanted to be able to experiment with an environment supporting functional data, such as mean functions in nonlinear regression models and prior density and likelihood functions in Bayesian analyses. Finally, I was interested in exploring the use of object oriented programming ideas for building and analyzing statistical models. I will discuss each of these points in a little more detail in the following paragraphs.
\piThe development of high resolution graphical computer displays has made it possible to consider the use of dynamic graphics for understanding higher-dimensional structure. One of the earliest examples is the real time rotation of a three dimensional point cloud on a screen -- an effort to use motion to recover a third dimension from a two dimensional display. Other techniques that have been developed include \Ibrushing\i a scatterplot -- highlighting points in one plot and seeing where the corresponding points fall in other plots. A considerable amount of research has been done in this area, see for example the discussion in Becker and Cleveland [4]. and the papers reproduced in Cleveland and McGill[8]. However most of the software developed to date has been developed on specialized hardware, such as the TTY 5620 terminal or Lisp machines. As a result, very few statisticians have had an opportunity to experiment with dynamic graphics first hand, and still fewer have had access to an environment that would allow them to implement dynamic graphics ideas of their own. Several commercial packages for microcomputers now contain some form of dynamic graphics, but most do not allow users to customize their plots or develop functions for producing specialized plots, such as dynamic residual plots. XLISP-STAT provides at least a partial solution to these problems. It allows the user to modify a scatter plot with Lisp level functions and provides means for modifying the way in which a plot responds to mouse actions. It is also easy to add functions written in C to the program. On the Macintosh this \pnhas to be done by adding to the source code. On some unix systems it is also possible to compile and dynamically load code written in C or FORTRAN.
\piAn integrated environment for statistical calculations and graphics is essential for developing an understanding of the uses of dynamic graphics in statistics and for developing new graphical techniques. Such an environment must essentially be a programming language. Its basic data types must include types that allow groups of numbers - data sets - to be manipulated as entire objects. But in model-based analyses numerical data are only part of the information being used. The remainder is the model itself. Sometimes a model is easily characterized by specifying a set of numbers. A normal linear regression model with \Ii. i. d.\i errors might be described by the number of covariates, the coefficients and the error variance. On the other hand, in many cases it is easier to specify a model by specifying a function. To specify a normal nonlinear regression model, for example, one might specify the mean function. If our language is to allow us to specify this function within the language itself then the language must support a functional data type with full rights: I\pnt has to be possible to define functions that manipulate functions, return functions, apply functions to arguments, etc.. The choice I faced was to define a language from scratch or use an existing language. Because of the complexity of issues involved in functional programming I decided to use a dialect of a well understood functional language, Lisp. The syntax of Lisp is somewhat unfamiliar to most users of statistical packages, but it is easy to learn and several good tutorials are available in local book stores. I considered the possibility of using Lisp to write a top level interface with a more ``natural'' syntax, but I did not see any way of doing this without complicating access to some of the more powerful features of Lisp or running into some of the pitfalls of functional programming. I therefore decided to retain the basic Lisp top level syntax. To make the manipulation of numerical data sets easier I have redefined the arithmetic operators and basic numerical functions to work on lists and arrays of data.
\piHaving decided to use Lisp as the basis for my environment XLISP was a natural choice for several reasons. It has been placed in the public domain by its author, David Betz. It is small (for a Lisp system), its source code is available in C, and it is easily extensible. Finally, it includes support for object oriented programming. Object oriented programming has received considerable attention in recent years and is particularly natural for use in describing and manipulating graphical objects. It may also be useful for the analysis of statistical data and models. A collection of data and assumptions may be represented as an \Iobject\i. The model object can then be examined and modified by sending it \Imessages\i. Many different kinds of models will answer similar questions, thus fitting naturally into an \Iinheritance structure\i. XLISP-STAT's implementation of linear and nonlinear regression models as \Iobjects\i, with nonlinear regression inheriting many of its\I methods\i from linear regression, is a first, primitive attempt to exploit this programming technique in statistical analysis.
\B\fs<14>Availability\b\fs<12>
Source code for XLISP-STAT for \IX11, Sun\i, 4.[23]BSD UNIX and Macintosh versions and executables for the Macintosh are available free of charge for non-commercial use. You should, however, be prepared to bear the cost of copying, for example by supplying a disk or tape and a stamped mailing envelope. You can also obtain the source code and Macintosh executables by anonymous \Iftp\i over the internet from \Iumnstat.stat.umn.edu\i (128.101.51.1).
\piA version for the Mac II requiring the MC68881 coprocessor is available as well. This version is compiled to access the coprocessor directly and will therefore not run on a machine without the coprocessor. Numerical computations with this version are about 7 to 10 times faster on a Mac II than without the direct coprocessor acces\pns.
\B\fs<14>Disclaimer\b\fs<12>
XLISP-STAT is an experimental program. It has not been extensively tested. The University of Minnesota, the School of Statistics, the author of the statistical extensions and the author of XLISP take no responsibility for losses or damages resulting directly or indirectly from the use of this program.
\piXLISP-STAT is an evolving system. Over time new features will be introduced, and existing features that do not work may be changed. Every effort will be made to keep XLISP-STAT consistent with the information in this tutorial, but if this is not possible the help information should give accurate information about the current use of a command.\pn
1\SpOn a UNIX system you can use the function \B\ff<(CG)Courier>load-data\b\ff<(CG)Times> to load the tutorial data. After starting up XLISP-STAT enter the expression\B\ff<(CG)Courier> (load-data "tutorial")\b\ff<(CG)Times>, followed by a \Ireturn\i.
\n\B\ff<(CG)Times>\fs<14>\jf\ll<1>\ps<100>\t<0>\Fp1 Starting and Finishing\b\fs<12>
You should have the program XLISP-STAT on a Macintosh disk. XLISP-STAT needs to have several files available for it to work properly. These files are\SP2\Sp:
\B\ff<(CG)Courier>init.lsp
common.lsp
help.lsp
objects.lsp
menus.lsp
statistics.lsp
dialogs.lsp
graphics.lsp
graphics2.lsp
regression.lsp
xlisp.help\b\ff<(CG)Times>\jf
Before starting XLISP-STAT you should make sure that these files are in the same folder as the XLISP-STAT application.
\piTo start XLISP-STAT double click on its icon. The program will need a little time to start up and read in the files mentioned above. When XLISP-STAT is ready the text in its command window will look something like this\SP3\Sp:
\B\ff<(CG)Courier>XLISP version 2.0, Copyright (c) 1988, by David Betz
XLISP-STAT version 2.0 , Copyright (c) 1988, by Luke Tierney.
Several files will be loaded; this may take a few minutes.
; loading "init.lsp"
; loading "common.lsp"
; loading "help.lsp"
; loading "objects.lsp"
; loading "menus.lsp"
; loading "statistics.lsp"
; loading "dialogs.lsp"
; loading "graphics.lsp"
; loading "graphics2.lsp"
; loading "regression.lsp"
>\b\ff<(CG)Times>
The final ">" in the window is the XLISP-STAT prompt. Any characters you type while the command window is active will be added to the line after the final prompt. When you hit a return, XLISP-STAT will try to interpret what you have typed and will print a response. For example, if you type a 1 and hit \Ireturn\i then XLISP-STAT will respond by simply printing a 1 on the following line and then giving you a new prompt:
\B\ff<(CG)Courier>
>\b\ff<(CG)Times>
\piIf you type an \Iexpression\i like \B\ff<(CG)Courier>(+ 1 2)\b\ff<(CG)Times>, then XLISP-STAT will print the result of evaluating the expression and give you a new prompt:\pn
\B\ff<(CG)Courier>
> (+ 1 2)
>\b\ff<(CG)Times>
As you have probably guessed, this expression means that the numbers 1 and 2 are to be added together. The next section will give more details on how XLISP-STAT expressions work. In this tutorial I will always show interactions with the program as I have done here: The ">" prompt will appear before lines you should type. XLISP-STAT will supply this prompt when it is ready; you should not type it yourself.
\piNow that you have seen how to start up XLISP-STAT it is a good idea to make sure you know how to get out. As with many Macintosh programs the easiest way to get out is to choose the \BQuit\b command from the \BFile\b menu. You can also use the command key shortcut COMMAND-Q, or you can type the expression\pn
\B\ff<(CG)Courier>
> (exit)\b\ff<(CG)Times>
Any one of these methods should cause the program to exit and return you to the Finder.
2\SpThe file \B\ff<(CG)Courier>xlisp.help\b\ff<(CG)Times> is optional. It may be replaced by a reduced file \B\ff<(CG)Courier>xlisp.help.small\b\ff<(CG)Times> or it may be omitted entirely. If it is not present interactive help will \pnnot be available.}:
\SP\pi3\SpOn a Macintosh with limited memory a dialog warning about memory restrictions may be appear at this point. On a Mac II it takes about a minute to load these files; on a Mac Plus or an SE it takes about 3.5 minutes.\fs<12>
\n\B\ff<(CG)Times>\fs<14>\jf\ll<1>\ps<100>\t<0>\Fp2 Introduction to Basics\b\fs<12>
Before we can start to use XLISP-STAT for statistical work we need to learn a little about the kind of data XLISP-STAT uses and about how the XLISP-STAT \Ilistener\i and \Ievaluator\i work.
\B2.1 Data\b
XLISP-STAT works with two kinds of data: \Isimple data\i and \Icompound data\i. Simple data are number\ll<-2>s
\B\ff<(CG)Courier>
1\s\ll<1>\s; an integer
-3.14\s\s; a floating point number
#C(0 \ll<-2>1)\s; a complex number (the imaginary unit)\b\ff<(CG)Times>
logical \ll<-5>v\ll<-2>alues
\B\ff<(CG)Courier>
T\s\ll<1>\s; true
nil\ll<-2>\s\s; false
\b\ff<(CG)Times>
string\ll<1>s (\ll<-2>always enclosed in double quotes)
\B\ff<(CG)Courier>"Thi\ll<-5>s\ll<-2> is a string 1 2 3 4"\b\ff<(CG)Times>
and \ll<-5>sy\ll<-2>mbols (used for naming things; see the following section)
\B\ff<(CG)Courier>
x\ll<1>
this\ll<-2>-is-a-symbol\b\ff<(CG)Times>
Compound data are lists
\B\ff<(CG)Courier>
(this \ll<1>is a list with 7 elements)
(+ 1 2 3)
(sq\ll<-2>rt 2)\b\ff<(CG)Times>
or \ll<-6>vecto\ll<-3>rs
\I\ff<(CG)Courier>
\B\i#(th\ll<1>is is a vector with 7 elements)
#(\ll<-2>1 2 3)\b\ff<(CG)Times>
Higher d\ll<1>imensional arrays are another form of compound data; they will be discussed below in Section 9.
\piAll the examples given above can be typed directly into the command window as they are shown here. The next section describes what XLISP-STAT will do with these expressions.\pn
\B2.2 The Listener and the Evaluator\b
A session with XLISP-STAT basically consists of a conversation between you and the \Ilistener\i. The listener is the window into which you type your commands. When it is ready to receive a command it gives you a prompt. At the prompt you can type in an expression. You can use the mouse or the \Ibackspace\i key to correct any mistakes you make while typing in your expression. When the expression is complete and you type a \Ireturn\i the listener passes the expression on to the \Ievaluator\i. The evaluator evaluates the expression and returns the result to the listener for printing.\SP4\Sp The evaluator it the heart of the system.
\piThe basic rule to remember in trying to understand how the evaluator works is that everything is evaluated. Numbers and strings evaluate to themselves:
\B\ff<(CG)Courier>\pn
> "Hello"
"Hello"
>\b\ff<(CG)Times>
Lists are more complicated. Suppose you type the list \B\ff<(CG)Courier>(+ 1 2 3)\b\ff<(CG)Times> at the listener. This list has four elements: the symbol \B\ff<(CG)Courier>+\b\ff<(CG)Times> followed by the numbers 1, 2 and 3. Here is what happens:
\B\ff<(CG)Courier>\ll<-2>
> (+\ll<1> 1 2 3)
>\b\ff<(CG)Times>\ll<-2>
\ll<1>
A list is evaluated as a function application. The first element is a symbol representing a function, in this case the symbol \B\ff<(CG)Courier>+\b\ff<(CG)Times> representing the addition function. The remaining elements are the arguments. Thus the list in the example above is interpreted to mean "Apply the function \B\ff<(CG)Courier>+\b\ff<(CG)Times> to the numbers 1, 2 and 3".
\piActually, the arguments to a function are always evaluated before the function is applied. In the previous example the arguments are all numbers and thus evaluate to themselves. On the other hand, consider\pn
\B\ff<(CG)Courier>\ll<-2>
> (\ll<1>+ (* 2 3) 4)
>\b\ff<(CG)Times>\ll<-2>
\ll<1>
The evaluator has to evaluate the first argument to the function \B\ff<(CG)Courier>+\b\ff<(CG)Times> before it can apply the function.
\piOccasionally you may want to tell the evaluator \Inot\i to evaluate something. For example, suppose we wanted to get the evaluator to simply return the list \B\ff<(CG)Courier>(+ 1 2)\b\ff<(CG)Times> back to us, instead of evaluating it. To do this we need to \Iquote\i our \pnlist:
\B\ff<(CG)Courier>\ll<-2>
> (q\ll<1>uote (+ 1 2))
(+ 1 2)
>\b\ff<(CG)Times>\ll<-2>
\ll<1>
\B\ff<(CG)Courier>quote\b\ff<(CG)Times> is not a function. It does not obey the rules of function evaluation described above: Its argument is not evaluated. \B\ff<(CG)Courier>quote\b\ff<(CG)Times> is called a \Ispecial form\i - special because it has special rules for the treatment of its arguments. There are a few other special forms that we will need; I will introduce them as they are needed. Together with the basic evaluation rules described here these special forms make up the basics of the LISP language. The special form \B\ff<(CG)Courier>quote\b\ff<(CG)Times> is used so often that a shorthand notation has been developed, a single quote before the expression you want to quote:
\B\ff<(CG)Courier>\ll<-3>
> '(+ 1 2) \s; single quote shorthand
\ll<1>
\b\ff<(CG)Times>This is equivalent to \B\ff<(CG)Courier>(quote (+ 1 2))\b\ff<(CG)Times>. Note that there is no matching quote following the expression.
\piBy the way, the semicolon ";" is the LISP comment character. Anything you type after a semicolon up to the next time you hit a return is ignored by the evaluator.\pn
\BExercises
For each of the following expressions try to predict what the evaluator will return. Then type them in, see what happens and try to explain any differences.
4\SpIt is possible to make a finer distinction. The \Ireader\i takes a string of characters from the listener and converts it into an expression. The \Ievaluator\i evaluates the expression and the \Iprinter\i converts the result into another string of characters for the listener to print. For simplicity I will use \Ievaluator\i to describe the combination of these functions.\fs<12>\pn
This section introduces some of the basic graphical and numerical statistical operations that are available in XLISP-STAT.
\B3.1 First Steps\b
Statistical data usually consists of groups of numbers. Devore and Peck [11, Exercise 2.11] describe an experiment in which 22 consumers reported the number of times they had purchased a product during the previous 48 week period. The results are given as a table:
\ll<-2>
\s\s\s\s0 \ll<1> 2 5 0 3 1 8 0 3 1 1
\s\s\s\s9\ll<-2> 2 4 0 2 9 3 0 1 9 8
\ll<1>
To examine this data in XLISP-STAT we represent it as a list of numbers using the \B\ff<(CG)Courier>list\b\ff<(CG)Times> function:
Note that the numbers are separated by white space (spaces, tabs or even returns), not commas.
\piThe \B\ff<(CG)Courier>mean\b\ff<(CG)Times> function can be used to compute the average of a list of numbers. We can combine it with the \B\ff<(CG)Courier>list\b\ff<(CG)Times> function to find the average number of purchases for our sample:\pn
It is of course a nuisance to have to type in the list of 22 numbers every time we want to compute a statistic for the sample. To avoid having to do this I will give this list a name using the \B\ff<(CG)Courier>def\b\ff<(CG)Times> special form\SP5\Sp:
Now the symbol \B\ff<(CG)Courier>purchases\b\ff<(CG)Times> has a value associated with it: Its value is our list of 22 numbers. If you give the symbol \B\ff<(CG)Courier>purchases\b\ff<(CG)Times> to the evaluator then it will find the value of this symbol and return that value:
\B\ff<(CG)Courier>\ll<-2>
> purch\ll<1>ases
(0 2 5 0 3 1 8 0 3 1 1 9 2 4 0 2 9 3 0 1 9 8)
> \b\ff<(CG)Times>
We can now easily compute various numerical descriptive statistics for this data set:
\B\ff<(CG)Courier>
> (mean purchases)
3.227273
> (median purchases)
> (standard-deviation purchases)
3.279544
> (interquartile-range purchases)
>\b\ff<(CG)Times>
\piXLISP-STAT also supports elementwise arithmetic operations on lists of numbers. For example, we can add 1 to each of the purchases:\pn
\B\ff<(CG)Courier>> (+ 1 purchases)
(1 3 6 1 4 2 9 1 4 2 2 10 3 5 1 3 10 4 1 2 10 9)
>\b\ff<(CG)Times>
and after adding 1 we can compute the natural logarithms of the results:
For each of the following expressions try to predict what the evaluator will return. Then type them in, see what happens and try to explain any differences.
Devore and Peck [11,page 54, Table 10] give precipitation levels recorded during the month of March in the Minneapolis - St. Paul area over a 30 year period. Let's enter these data into XLISP-STAT with the name \B\ff<(CG)Courier>precipitation\b\ff<(CG)Times>:
\jcFigure 1: Histogram of precipitation levels.\jf
In typing this expression I have hit the \Ireturn\i and \Itab\i keys a few times in order to make the typed expression easier to read. The tab key indents the next line to a reasonable point to make the
expression more readable.
\piThe \B\ff<(CG)Courier>histogram\b\ff<(CG)Times> and \B\ff<(CG)Courier>boxplot\b\ff<(CG)Times> functions can be used to obtain graphical representations of this data set:\pn
\B\ff<(CG)Courier>\ll<-2>
> (histogr\ll<1>am precipitation)
#<Object: 3564170, prototype = HISTOGRAM-PROTO>
> (boxplot precipitation)
#<Object: 3423466, prototype = SCATTERPLOT-PROTO>
>\b\ff<(CG)Times>\ll<-2>
\ll<1>
Each of these commands should cause a window with the appropriate graph to appear on your screen. The windows should look something Figures 1 and 2.
\piNote that as each graph appears it becomes the active window. To get XLISP-STAT to accept further commands you have to click on the XLISP-STAT listener window. You will have to click on the listener window between the two commands shown here.
The two functions return results that are printed something like this:\pn
These result will be used later to identify the window containing the plot. For the moment you can ignore them.
\piWhen you have several plot windows open you might want to close the listener window so you can rearrange the plots more easily. You can do this by clicking in the listener window's close box. You can later reopen the listener window by selecting the \BShow XLISP-STAT\b item on the \BCommand\b menu.\pn
\piHere are some numerical summaries:\pn
\B\ff<(CG)Courier>\ll<-2>
> (mea\ll<1>n precipitation)
1.685
> (median precipitation)
> (standard-deviation precipitation)
1.0157
> (interquartile-range precipitation)
1.145
\b\ff<(CG)Times>\jcFigure 2: Boxplot of precipitation levels.\jf
\piThe distribution of this data set is somewhat skewed to the right. Notice the separation between the mean and the median. You might want to try a few simple transformations to see if you can symmetrize the data. Square root and log transformations can be computed using the expressions \pn
You should look at plots of the data to see if these transformations do indeed lead to a more symmetric shape. The means and medians of the transformed data are
\B\ff<(CG)Courier>\ll<-2>
> (m\ll<1>ean (sqrt precipitation))
1.243006
> (median (sqrt precipitation))
1.212323
> (mean (log precipitation))
0.3405517
> (median (log precipitation))
0.384892
> \b\ff<(CG)Times>\ll<-2>
\ll<1>
The boxplot function can also be used to produce parallel boxplots of two or more samples. It will do so if it is given a list of lists as its argument instead of a single list. As an example, let's use this function to compare serum total cholesterol values for samples of rural and urban Guatemalans (see
Devore and Peck [11, page 19, Example 3]:
\jcFigure 3: Parallel box plots of cholesterol levels for urban and rural guatemalans.\jf
and is shown in Figure 3; the urban group is on the left.
\BExercises\b
The following exercises involve examples and problems from Devore and Peck [11]. The data sets are in files in the folder \BData\b on the XLISP-STAT distribution disk and can be read in using the \BLoad\b item in the \BFile\b menu or using the \B\ff<(CG)Courier>load\b\ff<(CG)Times> function (see Section 5.4 below). To use the \BLoad\b item on the \BFile\b menu select this item from the menu. This will bring up an \BOpen File\b dialog window. Use this dialog to open the \BData\b folder on the distribution disk. Now select one of the \B\ff<(CG)Courier>.lsp\b\ff<(CG)Times> files \B\ff<(CG)Courier>(car-prices.lsp\b\ff<(CG)Times> for the first exercise) and press the \BOpen\b button. The file will be loaded and some variables will be defined for you. Loading file \B\ff<(CG)Courier>car-prices.lsp\b\ff<(CG)Times> will define the single variable \B\ff<(CG)Courier>car-prices\b\ff<(CG)Times>. Loading file \B\ff<(CG)Courier>heating.lsp\b\ff<(CG)Times> will define two variables, \B\ff<(CG)Courier>gas-heat\b\ff<(CG)Times> and \B\ff<(CG)Courier>electric-heat\b\ff<(CG)Times>.\SP6\Sp
\po1. Devore and Peck [11, page 18, Example 2] give advertised prices for a sample of 50 used Japanese subcompact cars. Obtain some plots and summary statistics for this data. Experiment with some transformations of the data as well. The data set is called car-prices} in the file \B\ff<(CG)Courier>car-prices.lsp\b\ff<(CG)Times>. The prices are given in units of $1000; thus the price 2.39 represents $2390. The data have been sorted by their leading digit.\pn
\po2. In Exercise 2.40 Devore and Peck [11] give heating costs for a sample of apartments heated by gas and a sample of apartments heated by electricity. Obtain plots and summary statistics for these samples separately and look at a parallel box plot for the two samples. These data sets are called \B\ff<(CG)Courier>gas-heat\b\ff<(CG)Times> and \B\ff<(CG)Courier>electric-heat\b\ff<(CG)Times> in the file\B\ff<(CG)Courier> heating.lsp\b\ff<(CG)Times>.\pn
\B3.3 Two Dimensional Plots\b
Many single samples are actually collected over time. The precipitation data set used above is an example of this kind of data. In some cases it is reasonable to assume that the observations are independent of one another, but in other cases it is not. One way to check the data for some form of serial correlation or trend is to plot the observations against time, or against the order in which they were obtained. I will use the \B\ff<(CG)Courier>plot-points\b\ff<(CG)Times> function to produce a scatterplot of the precipitation data versus time. The \B\ff<(CG)Courier>plot-points\b\ff<(CG)Times> function is called as
\B\ff<(CG)Courier>\ll<-5>
(plot-points x-variable y-variable)
\b\ff<(CG)Times>\ll<1>
Our \Iy\i-variable will be \B\ff<(CG)Courier>precipitation\b\ff<(CG)Times>, the variable we defined earlier. As our \Ix\i-variable we would like to use a sequence of integers from 1 to 30. We could type these in ourselves, but there is an easier way. The function \B\ff<(CG)Courier>iseq\b\ff<(CG)Times>, short for \Iinteger-sequence\i, generates a list of consecutive integers between two specified values. The general form for a call to this function is
and the result will look like Figure 4. There does not appear to be much of a pattern to the data; an independence assumption may be reasonable.
\piSometimes it is easier to see temporal patterns in a plot if the points are connected by lines. Try the above command with \B\ff<(CG)Courier>plot-points\b\ff<(CG)Times> replaced by \B\ff<(CG)Courier>plot-lines\b\ff<(CG)Times>.\pn
\piThe \B\ff<(CG)Courier>plot-lines\b\ff<(CG)Times> function can also be used to construct graphs of functions. Suppose you would like a plot of sin(\Ix\i) from -\ff<(CG)Symbol>p\ff<(CG)Times> to +\ff<(CG)Symbol>p\ff<(CG)Times>. The constant \ff<(CG)Symbol>p\ff<(CG)Times> is predefined as the variable \B\ff<(CG)Courier>pi\b\ff<(CG)Times>. You can construct a list of\I n\i equally spaced real numbers between \Ia\i and \Ib\i using the expression\pn
\ll<-2>
\B\ff<(CG)Courier>(rseq a b n)\b\ff<(CG)Times>.
\ll<1>
Thus to draw the plot of sin(\Ix\i) using 50 equally spaced points type
\jcFigure 4: Scatterplot of precipitation levels against time.\jf
\B\ff<(CG)Courier>\ll<-2>
> (plot-lines\ll<1> (rseq (- pi) pi 50) (sin (rseq (- pi) pi 50)))
#<Object: 3423466, prototype = SCATTERPLOT-PROTO>
> \b\ff<(CG)Times>\ll<-2>
\ll<1>
The plot should look like Figure 5.
\piScatterplots are of course particularly useful for examining the relationship between two numerical observations taken on the same subject. Devore and Peck [11, Exercise 2.33] give data for HC and CO emission recorded for 46 automobiles. The results can be placed in two variables, \B\ff<(CG)Courier>hc\b\ff<(CG)Times> and \B\ff<(CG)Courier>co\b\ff<(CG)Times>, and these variable can then be plotted against one another with the \B\ff<(CG)Courier>plot-points\b\ff<(CG)Times> function:\pn
5\Spdef\b\ff<(CG)Times> acts like a special form, rather than a function, since its first argument is not evaluated (otherwise you would have to quote the symbol). Technically \B\ff<(CG)Courier>def\b\ff<(CG)Times> is a macro, not a special form, but I will not worry about this distinction in this document. \B\ff<(CG)Courier>def\b\ff<(CG)Times> is closely related to the standard Lisp special forms \B\ff<(CG)Courier>setf\b\ff<(CG)Times> and \B\ff<(CG)Courier>setq\b\ff<(CG)Times>. The advantage of using \B\ff<(CG)Courier>def\b\ff<(CG)Times> is that it adds your variable name\pn to a list of \B\ff<(CG)Courier>def\b\ff<(CG)Times>'ed variables kept in the global variable \B\ff<(CG)Courier>variables\b\ff<(CG)Times>. If you use \B\ff<(CG)Courier>setf\b\ff<(CG)Times> or \B\ff<(CG)Courier>setq\b\ff<(CG)Times> there is no easy way to find variables you have defined, as opposed to ones that are predefined. \B\ff<(CG)Courier>def\b\ff<(CG)Times> always affects top level symbol bindings, not local bindings. It can not be used in function definitions to change local bindings.\fs<12>
7\SpOn UNIX systems use the function \B\ff<(CG)Courier>load-data\b\ff<(CG)Times>. For example, evaluating the expression \B\ff<(CG)Courier>(load-data "car-prices")\b\ff<(CG)Times> should load the file\B\ff<(CG)Courier> car-prices.lsp\b\ff<(CG)Times>.
1. Draw a graph of the function f(x) = 2x + x\SP2\Sp between -2 and 3.
2. \poDevore and Peck [11, Exercise 4.2] give the age and CPK concentration, a measure of metabolic activity, recorded for 18 cross country skiers during a relay race. These data are in the variables \B\ff<(CG)Courier>age\b\ff<(CG)Times> and \B\ff<(CG)Courier>cpk\b\ff<(CG)Times> in the file \B\ff<(CG)Courier>metabolism.lsp\b\ff<(CG)Times>. Plot the data and describe any relationship you observe between age and CPK concentration.\pn
\B3.4 Plotting Functions\b
Plotting the sine functions in the previous section was a bit cumbersome. As an alternative we can use the function \B\ff<(CG)Courier>plot-function\b\ff<(CG)Times> to plot a function of one argument over a specified range. We can plot the sine function using the expression
\B\ff<(CG)Courier>
(plot-function (function sin) (- pi) pi)
\b\ff<(CG)Times>
The expression \B\ff<(CG)Courier>(function sin)\b\ff<(CG)Times> is needed to extract the function associated with the symbol \B\ff<(CG)Courier>sin\b\ff<(CG)Times>. Just using \B\ff<(CG)Courier>sin\b\ff<(CG)Times> will not work. The reason is that a symbol in Lisp can have both a \Ivalue\i, perhaps set using \B\ff<(CG)Courier>def\b\ff<(CG)Times>, and a \Ifunction definition\i at the same time.\SP7\Sp This may seem a bit cumbersome at first, but it has one great advantage: Typing an innocent expression like
\B\ff<(CG)Courier>
(def list '(2 3 4))\b\ff<(CG)Times>
will not destroy the\B\ff<(CG)Courier> list\b\ff<(CG)Times> function.
Extracting a function definition from a symbol is done almost as often as quoting an expression, so again a simple shorthand notation is available. The expression
\B\ff<(CG)Courier>#'sin
\b\ff<(CG)Times>is equivalent to the expression \B\ff<(CG)Courier>(function sin).\b\ff<(CG)Times>The short form \B\ff<(CG)Courier>#'\b\ff<(CG)Times> is usually pronounc-ed \Isharp-quote\i. Using this abbreviation the expression for producing the sine plot can be written as
\B\ff<(CG)Courier>
(plot-function #'sin (- pi) pi)
a\n\ff<(CG)Times>\fs<12>\jf\ll<1>\ps<100>\t<0>\Fp \SP\fs<10>\pi7\SpAs an aside, a Lisp symbol can be thought of as a "thing" with four cells. These cells contain the symbol's print name, its value, its function definition, and its property list. Lisp symbols are thus much more like physical entities than variable identifiers in FORTRAN or C.\pn
will generate a list of 50 independent uniform random variables. The functions \B\ff<(CG)Courier>normal-rand\b\ff<(CG)Times> and \B\ff<(CG)Courier>cauchy-rand\b\ff<(CG)Times> work similarly. Other generating functions require additional arguments to specify distribution parameters. Here is a list of the functions available for dealing with probability distributions:\ll<-2>
More informatio\ll<1>n on the required arguments is given in the appendix in Section C.3. The discrete quantile functions \B\ff<(CG)Courier>binomial-quant\b\ff<(CG)Times> and \B\ff<(CG)Courier>poisson-quant\b\ff<(CG)Times> return values of a left continuous inverse of the cdf. The pmf's for these distributions are only defined for integer arguments. The quantile functions and random variable generators for the beta, gamma, \ff<(CG)Symbol>c\SP\ff<(CG)Times>2\Sp, t and F distributions are presently calculated by inverting the cdf and may be a bit slow.
\piThe state of the internal random number generator can be ``randomly'' reseeded, and the current value of the generator state can be saved. The mechanism used is the standard Common Lisp mechanism. The current random state is held in the variable \B*\ff<(CG)Courier>random-state*\b\ff<(CG)Times>. The function \B\ff<(CG)Courier>make-random-state\b\ff<(CG)Times> can be used to set and save the state. It takes an optional argument. If the argument is \B\ff<(CG)Courier>NIL\b\ff<(CG)Times> or omitted \B\ff<(CG)Courier>make-random-state\b\ff<(CG)Times> returns a copy of the current value of \B\ff<(CG)Courier>*random-state*\b\ff<(CG)Times>. If the argument is a state object a copy of it is returned. If the argument is\B\ff<(CG)Courier> t\b\ff<(CG)Times> a new, "randomly" initialized state object is produced and returned.\SP8\Sp\pn
\B4.2 Generating Systematic Data\b
We have already used the functions \B\ff<(CG)Courier>iseq\b\ff<(CG)Times> and \B\ff<(CG)Courier>rseq\b\ff<(CG)Times> to generate equally spaced sequences of integers and real numbers. The function \B\ff<(CG)Courier>repeat\b\ff<(CG)Times> is useful for generating sequences with a particular pattern. The general form of a call to \B\ff<(CG)Courier>repeat\b\ff<(CG)Times> is
\B\ff<(CG)Courier>(repeat list pattern)\b\ff<(CG)Times>
\B\ff<(CG)Courier>pattern\b\ff<(CG)Times> must be either a single number or a list of numbers of the same length as \B\ff<(CG)Courier>list\b\ff<(CG)Times>. If \B\ff<(CG)Courier>pattern\b\ff<(CG)Times> is a single number then \B\ff<(CG)Courier>repeat\b\ff<(CG)Times> simply repeats \B\ff<(CG)Courier>list pattern\b\ff<(CG)Times> times. For example
\B\ff<(CG)Courier>> (repeat (list 1 2 3) 2)
(1 2 3 1 2 3)\b\ff<(CG)Times>\ll<-2>
\ll<1>
If \B\ff<(CG)Courier>pattern\b\ff<(CG)Times> is a list then each element of \B\ff<(CG)Courier>list\b\ff<(CG)Times> is repeated the number of times indicated by the corresponding element of \B\ff<(CG)Courier>pattern\b\ff<(CG)Times>. For example
In Section 6.2 below I generate the variables \B\ff<(CG)Courier>density\b\ff<(CG)Times> and \B\ff<(CG)Courier>variety\b\ff<(CG)Times> by typing them in directly. Using the \B\ff<(CG)Courier>repeat\b\ff<(CG)Times> function we could have generated them like this:
The \B\ff<(CG)Courier>select\b\ff<(CG)Times> function allows you to select a single element or a group of elements from a list or vector. For example, if we define \B\ff<(CG)Courier>x\b\ff<(CG)Times> by
then \B\ff<(CG)Courier>(select x i)\b\ff<(CG)Times> will return the\I i\i-th element of \B\ff<(CG)Courier>x\b\ff<(CG)Times>. LISP, like the language C but in contrast to Fortran, numbers elements of list and vectors starting at zero. Thus the indices for the elements of \B\ff<(CG)Courier>x\b\ff<(CG)Times> are 0, 1, 2, 3, 4, 5, 6, 7. So
\ll<-1>
\B\ff<(CG)Courier>> (select x\ll<1> 0)
> (select x 2)
5\b\ff<(CG)Times>\ll<-5>
\ll<1>
To get a group of elements at once we can use a list of indices instead of a single index:
\ll<-2>
\B\ff<(CG)Courier>> (sele\ll<1>ct x (list 0 2))
(3 5)\b\ff<(CG)Times>\ll<-2>
\ll<1>
\piIf you want to select all elements of x except element 2 you can use the expression \pn
Another approach is to use the logical function /= (meaning not equal) to tell you which indices are not equal to 2. The function \B\ff<(CG)Courier>which\b\ff<(CG)Times> can then be used to return a list of all the indices for which the elements of its argument are not \B\ff<(CG)Courier>NIL\b\ff<(CG)Times>:
\ll<-2>
\B\ff<(CG)Courier>> (/= 2 (iseq 0 7))\ll<1>
(T T NIL T T T T T)
> (which (/= 2 (iseq 0 7)))
(0 1 3 4 5 6 7)
> (select x (which (/= 2 (iseq 0 7))))
(3 7 9 12 3 14 2)\b\ff<(CG)Times>
This approach is a little more cumbersome for deleting a single element, but it is more general. The expression \B\ff<(CG)Courier>(select x (which (< 3 x)))\b\ff<(CG)Times>, for example, returns all elements in \B\ff<(CG)Courier>x\b\ff<(CG)Times> that are greater than 3:\ll<-2>
At times you may want to combine several short lists into a single longer list. This can be done using the \B\ff<(CG)Courier>append\b\ff<(CG)Times> function. For example, if you have three variables \B\ff<(CG)Courier>x\b\ff<(CG)Times>, \B\ff<(CG)Courier>y\b\ff<(CG)Times> and \B\ff<(CG)Courier>z\b\ff<(CG)Times> constructed by the expressions\ll<-2>
So far when I have asked you to type in a list of numbers I have been assuming that you will type the list correctly. If you made an error you had to retype the entire\B\ff<(CG)Courier> def\b\ff<(CG)Times> expression. Since you can use cut--and--paste this is really not too serious. However it would be nice to be able to replace the values in a list after you have typed it in. The \B\ff<(CG)Courier>setf\b\ff<(CG)Times> special form is used for this. Suppose you would like to chang\ll<-2>e the 12 in the list \B\ff<(CG)Courier>x\b\ff<(CG)Times> used in the Section 4.3 to 11. The expression
\ll<1>
\B\ff<(CG)Courier>(setf (se\t<1>\ll<-2>lect x 4) 11)\b\ff<(CG)Times>
will m\ll<-5>ake thi\t<0>\ll<-2>s replacement:
\B\ff<(CG)Courier>> (s\ll<1>etf (select x 4) 11)
(3 7 5 \ll<-2>9 11 3 14 2)\b\ff<(CG)Times>
\piTh\ll<1>e general form of \B\ff<(CG)Courier>setf\b\ff<(CG)Times> is\pn
\ll<-2>
\B\ff<(CG)Courier>(setf form value)\b\ff<(CG)Times>
where \B\ff<(CG)Courier>form\b\ff<(CG)Times> i\ll<1>s the expression you would use to select a single element or a group of elements from \B\ff<(CG)Courier>x\b\ff<(CG)Times> and \B\ff<(CG)Courier>value\b\ff<(CG)Times> is the value you would like that element to have, or the list of the values for the elements in the group. Thus the expression
\ll<-2>
\B\ff<(CG)Courier>(setf (select x (list 0 2)) (list 15 16))\b\ff<(CG)Times>
\ll<1>
changes the values \ll<-2>of elements 0 and 2 to 15 and 16:
\B\ff<(CG)Courier>> (setf (sele\ll<1>ct x (list 0 2)) (list 15 16))
(15 16)
(15 7 16 9 11 3 14 2)\b\ff<(CG)Times>
\piA note of caution is needed here. Lisp symbols are merely labels for different items. When you assign a name to an item with the \B\ff<(CG)Courier>def\b\ff<(CG)Times> command you are not producing a new item. Thus
\B\ff<(CG)Courier>(def x (list 1 2 3 4))
(def y x)\b\ff<(CG)Times>
means that \B\ff<(CG)Courier>x\b\ff<(CG)Times> and \B\ff<(CG)Courier>y\b\ff<(CG)Times> are two different names for the same thing. As a result, if we change an element of (the item referred to by) \B\ff<(CG)Courier>x\b\ff<(CG)Times> with \B\ff<(CG)Courier>setf\b\ff<(CG)Times> then we are also changing the element of (the item referred to by) \B\ff<(CG)Courier>y\b\ff<(CG)Times>, since both \B\ff<(CG)Courier>x\b\ff<(CG)Times> and\B\ff<(CG)Courier> y\b\ff<(CG)Times> refer to the same item. If you want to make a copy of\B\ff<(CG)Courier> x\b\ff<(CG)Times> and store it in \B\ff<(CG)Courier>y\b\ff<(CG)Times> before you make changes to \B\ff<(CG)Courier>x\b\ff<(CG)Times> then you must do so explicitly using, say, the \B\ff<(CG)Courier>copy-list\b\ff<(CG)Times> function. The expression
\B\ff<(CG)Courier>(def y (copy-list x))\b\ff<(CG)Times>
will make a copy of \B\ff<(CG)Courier>x\b\ff<(CG)Times> and set the value of\B \ff<(CG)Courier>y\b\ff<(CG)Times> to that copy. Now \B\ff<(CG)Courier>x\b\ff<(CG)Times> and \B\ff<(CG)Courier>y\b\ff<(CG)Times> refer to different objects and changes to \B\ff<(CG)Courier>x\b\ff<(CG)Times> will not affect \B\ff<(CG)Courier>y\b\ff<(CG)Times>.
8\SpThe generator used is Marsaglia's portable generator from the \ICore Math Libraries\i distributed by the National Bureau of Standards. A state object is a vector containing the state information of the generator. "Random" reseeding occurs off the system clock.
Jt\n\B\ff<(CG)Times>\fs<14>\jf\ll<1>\ps<100>\t<0>\Fp5 Some Useful Shortcuts\b\fs<12>
This section describes some additional features of XLISP-STAT that you may find useful.
\B5.1 Getting Help\b
On line help is available for many of the functions in XLISP-STAT\SP9\Sp As an example, here is how you would get help for the function \B\ff<(CG)Courier>median\b\ff<(CG)Times>:
\B\ff<(CG)Courier>> (help 'median)
MEDIAN\s\s\s\s\s\s\s\s\s[function-doc]
Args: (x)
Returns the median of the elements of X.
> \b\ff<(CG)Times>
Note the quote in front of \B\ff<(CG)Courier>median\b\ff<(CG)Times>. \B\ff<(CG)Courier>help\b\ff<(CG)Times> is itself a function, and its argument is the symbol representing the function \B\ff<(CG)Courier>median\b\ff<(CG)Times>. To make sure \B\ff<(CG)Courier>help\b\ff<(CG)Times> receives the symbol, not the value of the symbol, you need to quote the symbol.
\piIf you are not sure about the name of a function you may still be able to get some help. Suppose you want to find out about functions related to the normal distribution. Most such functions will have "norm" as part of their name. The expression\pn
\B\ff<(CG)Courier>(help* 'norm)\b\ff<(CG)Times>
will print the help information for all symbols whose names contain the string "norm":\ff<(CG)Courier>>
The s\ll<1>ymbols \B\ff<(CG)Courier>norm\b\ff<(CG)Times>, \B\ff<(CG)Courier>norm-2\b\ff<(CG)Times> and \B\ff<(CG)Courier>normal\b\ff<(CG)Times> do not have function definitions and thus there is no help available for them. The term \Ivectorized\i in a function's documentation means the function can be applied to arguments that are lists or arrays; the result will be a list or array of the results of applying the function to each element of its arguments.\SP10\Sp A related term appearing in the documentation of some functions is \Ivector reducing\i. A function is vector reducing if it is applied recursively to its arguments until a single number results. The functions \B\ff<(CG)Courier>sum\b\ff<(CG)Times>, \B\ff<(CG)Courier>prod\b\ff<(CG)Times>, \B\ff<(CG)Courier>max\b\ff<(CG)Times> and \B\ff<(CG)Courier>min\b\ff<(CG)Times> are vector reducing.
\piThe first time a help function is used will take some time - the help file is scanned and the positions of all entries in the file are recorded. Subsequent calls will be faster.
Let me briefly explain the notation used in the information printed by \B\ff<(CG)Courier>help\b\ff<(CG)Times> to describe the arguments a function expects\SP11\Sp. Most functions expect a fixed set of arguments, described in the help message by a line like\pn\ll<-2>
\B\ff<(CG)Courier>Args: (x y z)\b\ff<(CG)Times>
Some functions can tak\ll<1>e one or more optional arguments. The arguments for such a function might be listed as\ll<-2>
\B\ff<(CG)Courier>Args: (x &optional y (z t))\b\ff<(CG)Times>
This means t\ll<1>hat \B\ff<(CG)Courier>x\b\ff<(CG)Times> is required and \B\ff<(CG)Courier>y\b\ff<(CG)Times> and \B\ff<(CG)Courier>z\b\ff<(CG)Times> are optional. If the function is named \B\ff<(CG)Courier>f\b\ff<(CG)Times>, it can be called as \B\ff<(CG)Courier>(f x-val)\b\ff<(CG)Times>, \B\ff<(CG)Courier>(f x-val y-val)\b\ff<(CG)Times> or \B\ff<(CG)Courier>(f x-val y-val z-val)\b\ff<(CG)Times>. The list \B\ff<(CG)Courier>(z t)\b\ff<(CG)Times> means that if \B\ff<(CG)Courier>z\b\ff<(CG)Times> is not supplied its default value is \B\ff<(CG)Courier>T\b\ff<(CG)Times>. No explicit default value is specified for \B\ff<(CG)Courier>y\b\ff<(CG)Times>; its default value is therefore \B\ff<(CG)Courier>NIL\b\ff<(CG)Times>. The arguments must be supplied in the order in which they are listed. Thus if you want to give the argument \B\ff<(CG)Courier>z\b\ff<(CG)Times> you must also give a value for \B\ff<(CG)Courier>y\b\ff<(CG)Times>.
\piAnother form of optional argument is the \Ikeyword argument\i. Thehistogram function for example takes arguments\pn\ll<-3>
The \B\ff<(CG)Courier>data\b\ff<(CG)Times> argument is req\ll<0>uir\ll<1>ed, the \B\ff<(CG)Courier>title\b\ff<(CG)Times> argument is an optional keyword argument. The default title is \B\ff<(CG)Courier>"Histogram"\b\ff<(CG)Times>. If you want to create a histogram of the purchases data set used in Section 3.1 with title \B\ff<(CG)Courier>"Purchases"\b\ff<(CG)Times> use the expression\ll<-3>
Thus to give a value for a keyword argument\ll<0> \ll<1>you give the keyword\SP12\Sp for the argument, a symbol consisting of a colon followed by the argument name, and then the value for the argument. If a function can take several keyword arguments then these may be specified in any order, following the required and optional arguments.
\piFinally, some functions can take an arbitrary number of arguments. This is denoted by a line like\pn\ll<-3>
Args: (x &rest args)
The argum\ll<1>ent \B\ff<(CG)Courier>x\b\ff<(CG)Times> is required, and zero or more additional arguments can be supplied.
\piIn addition to providing information about functions \B\ff<(CG)Courier>help\b\ff<(CG)Times> also gives information about data types and certain variables. For example,\pn
\B\ff<(CG)Courier>> (help 'complex)
COMPLEX\s\s \s\s\s\s\s\s\s[function-doc]
Args: (realpart &optional (imagpart 0))
Returns a complex number with the given real and imaginary parts.
COMPLEX\s\s \s\s\s\s\s\s\s[type-doc]
A complex number
>\b\ff<(CG)Times>\ll<-3>
\ll<-2>
shows the function and type documentation for \B\ff<(CG)Courier>complex\b\ff<(CG)Times>, and
\ll<-3>
\B\ff<(CG)Courier>> (help\ll<0> \ll<1>'pi)
PI\s\s [variable-doc]
T\t<-1>he floating-point number that is approximately equal to the ratio of the\t<0>
circumference of a circle to its diameter.
>\b\ff<(CG)Times>\ll<-3>
\ll<1>
shows the variable documentation for \B\ff<(CG)Courier>pi\b\SP\ff<(CG)Times>13\Sp.
\B5.2 Listing and Undefining Variables\b
After you have been working for a while you may want to find out what variables you have defined (using \B\ff<(CG)Courier>def\b\ff<(CG)Times>). The function \B\ff<(CG)Courier>variables\b\ff<(CG)Times> will produce a listing:
\ll<-3>
\B\ff<(CG)Courier>> (varia\ll<0>b\ll<1>les)
RURAL
URBAN
PRECIPITATION
PURCHASES
>\b\ff<(CG)Times>\ll<-4>
\ll<-1>
\pi\ll<1>If you are working on a 1Mb Macintosh you may occasionally want to free up some space by getting rid of some variables you no longer need. You can do this using the\B \ff<(CG)Courier>undef\b\ff<(CG)Times> function:\pn
\ll<-3>
\B\ff<(CG)Courier>> (\ll<0>unde\ll<1>f 'co)
> (variables)
RURAL
URBAN
PRECIPITATION
PURCHASES
>\b\ff<(CG)Times>\ll<-3>
\ll<1>
\B5.3 More on the XLISP-STAT Listener\b\ll<0>
\ll<1>
Because of the large number of parentheses involved, Lisp expressions can be hard to read and type correctly. To make it easier to type readable, correct expressions the listener window on the Macintosh has the following features:
\s\po
Typing a left parenthesis flashes the matching right parenthesis.\ll<-1>
\ll<1>
You can move the cursor with the arrow keys, the mouse or the \Ibackspace\i key to any \T\Tposition in the current input expression, not just within the last line.\ll<-1>
\ll<1>
If the current expression is more than one line long, hitting the \Itab\i key in any line but the \T\Tfirst indents the line to its appropriate position according to (more or less) standard rules for \T\TLisp code indentation.\ll<-1>
\ll<1>
If the insertion point is at the end of the current expression hitting the \Ienter\i key is \T\Tequivalent to hitting a \Ireturn\i. If the insertion point is not at the end of the expression, \T\Thitting \Ienter\i moves the insertion point to the end of the expression.\pn
These four features should make typing expressions correctly much easier. In particular, in translating mathematical formulas to Lisp it sometimes seems that you have to do things backwards. Using these features you can build up your expression from the inside out\SP14\Sp.
\piXLISP 2.0 provides a simple command history mechanism. The symbols \B\ff<(CG)Courier>-\b\ff<(CG)Times>, \B\ff<(CG)Courier>*\b\ff<(CG)Times>, \B\ff<(CG)Courier>**\b\ff<(CG)Times>, \B\ff<(CG)Courier>***\b\ff<(CG)Times>, \B\ff<(CG)Courier>+\b\ff<(CG)Times>, \B\ff<(CG)Courier>++\b\ff<(CG)Times>, and \B\ff<(CG)Courier>+++\b\ff<(CG)Times> are used for this purpose. The top level reader binds these symbols as follows:\pn\ll<-2>
\ll<1>
\s\s\s\s\s \B\ff<(CG)Courier>-\b\ff<(CG)Times>\sthe current input expression
\s\s\s\s\s \B\ff<(CG)Courier>+\b\ff<(CG)Times>\sthe last expression read
\s\s\s\s\s \B\ff<(CG)Courier>++\s\b\ff<(CG)Times>the previous value of \B\ff<(CG)Courier>+\b\ff<(CG)Times>
\s\s\s\s\s\B\ff<(CG)Courier>+++\t<2>\s\b\ff<(CG)Times>\t<0>the previous value of \B\ff<(CG)Courier>++\b\ff<(CG)Times>
\s\s\s\s\s \B\ff<(CG)Courier>*\t<-1>\s\b\ff<(CG)Times>\t<0>the result of the last evaluation
\s\s\s\s\s \B\ff<(CG)Courier>**\s\b\ff<(CG)Times>the previous value of \B\ff<(CG)Courier>*\b\ff<(CG)Times>
\s\s\s\s\s\B\ff<(CG)Courier>***\s\b\ff<(CG)Times>the previous value of \B\ff<(CG)Courier>**\b\ff<(CG)Times>\ll<-2>
\ll<1>
The variables \B\ff<(CG)Courier>*\b\ff<(CG)Times>, \B\ff<(CG)Courier>**\b\ff<(CG)Times> and \B\ff<(CG)Courier>***\b\ff<(CG)Times> are probably most useful. For example, if you construct a plot but forget to assign the resulting plot object to a variable you can recover it using one of the history variables:\ll<-4>
\ll<1>
\B\ff<(CG)Courier>> (histogram (normal-rand 50))
#<Object: 3701682, prototype = HISTOGRAM-PROTO>
> (def w *)
#<Object: 3701682, prototype = HISTOGRAM-PROTO>
> \b\ff<(CG)Times>\ll<-4>
\ll<1>
The symbol \B\ff<(CG)Courier>W\b\ff<(CG)Times> now has the histogram object as its value and can be used to modify the plot, as described in Section 6.5 below.
\piLike most interactive systems, XLISP needs a system for dynamically managing memory. The system used by XLISP is to grab memory out of a fixed bin until the bin is exhausted. At that point the system pauses to reclaim memory that is no longer being used. This process, called \Igarbage collection\i, will occasionally cause the system to freeze up for about a second. When the system garbage collects the Macintosh \pncursor changes to a trash bag.
\piOccasionally a calculation will take too long, or it will appear to \pnhave gotten stuck in some kind of loop. If you want to interrupt the calculation \Ihold down\i the COMMAND key and the PERIOD. This should return you to the listener. You must continue to hold down the key until the calculation stops.
\B5.4 Loading Files\b
The data for the examples and exercises in this tutorial have been stored on files with names ending in \B\ff<(CG)Courier>.lsp\b\ff<(CG)Times>. On the XLISP-STAT distribution disk they can be found in the folder \B\ff<(CG)Courier>Data\b\ff<(CG)Times>. Any variables you save (see the next subsection for details) will also be saved on files of this form. The data in these files can be read into XLISP-STAT with the \B\ff<(CG)Courier>load\b\ff<(CG)Times> function. To load a file named \B\ff<(CG)Courier>randu.lsp\b\ff<(CG)Times> type the expression
If you give \B\ff<(CG)Courier>load\b\ff<(CG)Times> a name that does not end in\B \ff<(CG)Courier>.lsp\b\ff<(CG)Times> then load will add this suffix. Alternatively, you can use the \BLoad\b command in the \BFile\b menu. After loading a file using the \BFile\b menu \BLoad\b item the system does not print a prompt. Instead it prints a message indicating that the load is done:
starts a recording. All expressions typed by you and all results typed by XLISP-STAT will be entered into the file named \B\ff<(CG)Courier>myfile\b\ff<(CG)Times>. The expression
\ll<-2>
\B\ff<(CG)Courier>(dribble)\b\ff<(CG)Times>
\ll<1>
stops the recording. Note that \B\ff<(CG)Courier>(dribble "myfile")\b\ff<(CG)Times> starts a new file by the name \B\ff<(CG)Courier>myfile\b\ff<(CG)Times>. If you already have a file by that name its contents will be lost. Thus you can't use dribble to toggle on and off recording to a single file. You can also turn dribbling on and off using the \BDribble\b item on the \BCommand\b menu.
\B\ff<(CG)Courier>\pidribble\b\ff<(CG)Times> only records text that is typed, not plots. However, you can use the standard Macintosh shortcut COMMAND-SHIFT-3 to save a MacPaint image of the current screen. You can also choose the \BCopy\b command from the \BEdit\b menu, or its command key shortcut COMMAND-C, while a plot window is the active window to save the contents of the plot window to the clip board. You can then open the scrap book from the apple menu and paste the plot into the scrap book.
Variables you define in XLISP-STAT only exist for the duration of the current session. If you quit from XLISP-STAT, or the program crashes, your data will be lost. To preserve your data you can use the \B\ff<(CG)Courier>savevar\b\ff<(CG)Times> function. This function allows you to save one or more variables into a file. Again a new file is created and any existing file by the same name is destroyed. To save the variable \B\ff<(CG)Courier>precipitation\b\ff<(CG)Times> in a file called \B\ff<(CG)Courier>precipitation.lsp\b\ff<(CG)Times> type\pn
\b\ff<(CG)Times>\ll<1>Do not add the \B\ff<(CG)Courier>.lsp\b\ff<(CG)Times> suffix yourself; \B\ff<(CG)Courier>savevar\b\ff<(CG)Times> will supply it. To save the two variables \B\ff<(CG)Courier>precipitation\b\ff<(CG)Times> and \B\ff<(CG)Courier>purchases\b\ff<(CG)Times> in the file \B\ff<(CG)Courier>examples.lsp\b\ff<(CG)Times> type \SP15\Sp.
The files \B\ff<(CG)Courier>precipitation.lsp\b\ff<(CG)Times> and \B\ff<(CG)Courier>examples.lsp\b\ff<(CG)Times> now contain a set of expression that, when read in with the \B\ff<(CG)Courier>load\b\ff<(CG)Times> command, will recreate the variables \B\ff<(CG)Courier>precipitation\b\ff<(CG)Times> and \B\ff<(CG)Courier>purchases\b\ff<(CG)Times>. You can look at these files with an editor like MacWrite or the XLISP-STAT editor and you can prepare files with your own data by following these examples.
\B5.6 The XLISP-STAT Editor\b
The Macintosh version of XLISP-STAT includes a simple editor for preparing data files and function definitions. To edit an existing file select the \BOpen Edit\b item on the \BFile\b menu; to start a new file select \BNew Edit\b. Other commands on the \BFile\b menu can be used to save your file and to revert back to the saved version. The editor can only handle text files of less than 32K characters. As in the listener window, hitting the \Itab\i key in any line but the first of a multiline expression will indent the expression to a reasonable point. The editor also allows you to select a section of text representing one or more Lisp expressions and have these evaluated. To do this select the expressions you want to
evaluate and then choose \BEval Selection\b from the \BEdit\b menu. The returned values are not available, so this is only useful for producing side effects, such as defining variables or functions.
\B5.7 Reading Data Files\b
The data files we have used so far in this tutorial have been processed to contain XLISP-STAT expressions. XLISP-STAT also provides two functions for reading raw data files. The simpler of the two is \B\ff<(CG)Courier>read-data-file\b\ff<(CG)Times>. The expression
where\B \ff<(CG)Courier>file\b\ff<(CG)Times> is a string representing the name of the data file, returns a list of all the items in the file. Items can be separated by any amount of white space, but not by commas or other punctuation marks. Items can be any valid Lisp expressions. In particular they can be numbers, strings or symbols. The list can then be manipulated into the appropriate form within XLISP-STAT.
\piThe function \B\ff<(CG)Courier>read-data-columns\b\ff<(CG)Times> is provided for reading data files in which each row represents a case and each column a variable. The expression\pn
will return a list of \B\ff<(CG)Courier>cols\b\ff<(CG)Times> lists, each representing a column of the file. Note that this function determines the column structure from the value of \B\ff<(CG)Courier>cols\b\ff<(CG)Times>, not from the structure of the file. The entries of \B\ff<(CG)Courier>file\b\ff<(CG)Times> can be as for\B \ff<(CG)Courier>read-data-file\b\ff<(CG)Times>.
\piThese two functions should be adequate for most purposes. If you have to read a file that does not fit into the form considered here you can use the raw file handling functions of XLISP.\pn
\B5.8 User Initialization File\b
After loading in all its program files and before giving you your first prompt XLISP-STAT looks to see if there is a file named \B\ff<(CG)Courier>statinit.lsp\b\ff<(CG)Times> in the startup folder. If there is one it will be loaded. You can use this file to load any data sets you would like to have available or to define functions of your own.
9\Sp The on line help file may not be available on a single disk version that includes a system file. Alternatively, there may be a small help file available that does not contain documentation for all functions..\fs<12>
\SP\fs<10>\pi10\SpThis process of applying a function elementwise to its arguments is called \Imapping\i.\pn
\SP\pi11\SpThe notation used corresponds to the specification of the argument lists in Lisp function definitions. See Section 8 for more information on defining functions.\fs<12>
\SP\fs<10>12\Sp\pnNote that the keyword \B\ff<(CG)Courier>:title\b\ff<(CG)Times> has not been quoted. \IKeyword symbols\i, symbols starting with a colon, are somewhat special. When a keyword symbol is created its value is set to itself. Thus a keyword symbol effectively evaluates to itself and does not need to be quoted.\fs<12>\pi
\SP\pi14\SpTo support these features the listener checks the current expression each time you type a {return} to see if it has a complete expression. If so, the expression is passed to the reader and the evaluator. If not, you can continue typing. There are some heuristics involved here, and an expression with lots of quotes and comments may cause trouble, but it seems to work. Redefining the read table in major\fs<12> \fs<10>ways may not work as you think since some knowledge of standard Lisp syntax is built in to the listener.\fs<12>\pn
\SP\pi15\SpI have used a quoted list \B\ff<(CG)Courier>'(purchases precipitation)\ff<(CG)Times> \bin this expression to pass the list of symbols to the save function. A longer alternative would be the expression \B\ff<(CG)Courier>(list 'purchases 'precipitation)\b\ff<(CG)Times>.\pn
j&\n\B\ff<(CG)Times>\fs<14>\jf\ll<1>\ps<100>\t<0>\Fp6 More Elaborate Plots\b\fs<12>
The plots described so far were designed for exploring the distribution of a single variable and the relationship among two ariables. This section describes some plots and techniques that are useful for exploring the relationship among three or more variables. The techniques and plots described in the first four subsections are simple, requiring only simple commands to the listener and mouse actions. The fifth subsection introduces the concept of a \Iplot object\i and the idea of \Isending messages\i to an object from the listener window. These ideas are used to add lines and curves to scatter plots, and the basic concepts of objects and messages will be used again in the next section on regression models. The final subsection shows how Lisp iteration can be used to construct a dynamic simulation--a plot movie.
\B6.1 Spinning Plots\b
If we are interested in exploring the relationship among three variables then it is natural to imagine constructing a three dimensional scatterplot of the data. Of course we can only see a two dimensional projection of the plot on a computer screen -- any depth that you might be able to perceive by looking at a wire model of the data is lost. One approach to try to recover some of this depth perception is to rotate the points around some axis. The \B\ff<(CG)Courier>spin-plot\b\ff<(CG)Times> function allows you to construct a rotatable three dimensional plot.
\piAs an example let's look a some data collected to examine the relationship between a phosphate absorption index and the amount of extractable iron and aluminum in a sediment (Devore and Peck [11, page 509, Example 6]). The data can be entered with the expressions \pn
\B\ff<(CG)Courier>(spin-plot (list absorption iron aluminum))\b\ff<(CG)Times>
\ll<1>
produces the plot on the left in Figure 7. The argument to \B\ff<(CG)Courier>spin-plot\b\ff<(CG)Times> is a list of three lists or vectors, representing the \Ix\i, \Iy\i and \Iz\i variables. The \B\Iz\b\i axis is initially pointing out of the screen. You can rotate the plot by pointing at one of the \B\ff<(CG)Courier>Pitch\b\ff<(CG)Times>,\B \ff<(CG)Courier>Roll\b\ff<(CG)Times> or \B\ff<(CG)Courier>Yaw\b\ff<(CG)Times> squares and pressing the mouse button. By rotating the plot you can see that the points seem to fall close to a plane. The plot on the right of Figure 7 shows the data viewed along the plane. A linear model should describe this data quite well.
\piAs a second example, with the data defined by \pn
produces a plot that can be rotated to produce the view in Figure 8. These data concern samples of thawed permafrost soil. \B\ff<(CG)Courier>strength\b\ff<(CG)Times> is the shear strength, and \B\ff<(CG)Courier>water\b\ff<(CG)Times> is the percentage water content. \B\ff<(CG)Courier>depth\b\ff<(CG)Times> is the depth at which the sample was taken. The plot shows that a linear model will not fit well. Devore and Peck [11] suggest fitting a model with quadratic terms to this data.
\piThe function \B\ff<(CG)Courier>spin-plot\b\ff<(CG)Times> also accepts the additional keyword argument \B\ff<(CG)Courier>scale\b\ff<(CG)Times>. If \B\ff<(CG)Courier>scale\b\ff<(CG)Times> is \B\ff<(CG)Courier>T\b\ff<(CG)Times>, the default, then the data are centered at the midranges of the three variables, and all three variables are scaled to fit the plot. If scale is \B\ff<(CG)Courier>NIL\b\ff<(CG)Times> the data are assumed to be scaled between -1 and 1, and the plot is rotated about the origin. Thus if you want to center your plot at the means of the variables and scale all observations by the same amount you can use the expression\pn
Note that the \B\ff<(CG)Courier>scale\b\ff<(CG)Times> keyword argument is given using the corresponding keyword symbol, the symbol \B\ff<(CG)Courier>scale\b\ff<(CG)Times> preceded by a colon\pi.
\T \pnRotation speed can be changed using the plot menu or the keyboard equivalents COMMAND-F for Faster and COMMAND-S for Slower.\pi
Depth cuing and showing of the axes are controlled by the items on the plot menu.
If you click the mouse in one of the \B\ff<(CG)Courier>Pitch\b\ff<(CG)Times>, \B\ff<(CG)Courier>Roll\b\ff<(CG)Times> or \B\ff<(CG)Courier>Yaw\b\ff<(CG)Times> squares while holding down the \Ishift\i key the plot will start to rorate and continue to rotate after the mouse button has been released.\pn
\BExercises\b\ll<-2>
\ll<1>
\po1. An Upstate New York business machine company used to include a random number generator called RANDU in its software library. This generator was supposed to produce numbers that behaved like \Ii. i. d.\i uniform random variables. The data set randu in the file\B \ff<(CG)Courier>randu.lsp\b\ff<(CG)Times> in the \BData\b folder consists of a list of three lists of numbers. These lists are consecutive triples of numbers produced by RANDU. Use \B\ff<(CG)Courier>spin-plot\b\ff<(CG)Times> to see if you can spot any unusual features in this sample.\pn
\B6.2 \b \BScatterplot Matrices\b\ll<-2>
\ll<1>
Another approach to graphing a set of variables is to look at a matrix of all possible pairwise scatterplots of the variables. The\B \ff<(CG)Courier>scatterplot-matrix\b\ff<(CG)Times> function will produce such a plot. The data
produces the scatterplot matrix in Figure 9. The plot of\B\ff<(CG)Courier> abrasion-loss\b\ff<(CG)Times> against t\B\ff<(CG)Courier>ensile-strength\ff<(CG)Times> \bgives you an idea of the joint variation in these two variables. But \B\ff<(CG)Courier>hardness\b\ff<(CG)Times> varies from point to point as well. To get an understanding of the relationship among all three variables it would be nice to be able to fix \ff<(CG)Courier>hardness\ff<(CG)Times> at various levels and look at the way the plot of \B\ff<(CG)Courier>abrasion-loss\b\ff<(CG)Times> against \B\ff<(CG)Courier>tensile-strength\b\ff<(CG)Times> changes as you change these levels. You can do this kind of exploration in the scatterplot matrix by using the two highlighting techniques \B\ff<(CG)Courier>selecting\b\ff<(CG)Times> and \B\ff<(CG)Courier>brushing\b\ff<(CG)Times>.
\jcFigure 9: Scatterplot matrix of abrasion loss data.\jf
\s\I\poSelecting\i. Your plot is in the selecting mode when the cursor is an arrow. This is the default setting. In this mode select a point by clicking the mouse on top of it. To select a group of points drag a selection rectangle around the group. If the group does not fit in a rectangle you can build up your selection by holding down the \Ishift\i key as you click or drag. If you click without the \Ishift \ikey any existing selection will be unselected; when the \Ishift\i key is down selected points remain selected.\pn
\s\I\poBrushing\i. You can enter the brushing mode by choosing \BMouse Mode...\b from the \BScatmat\b menu and selecting \BBrushing\b from the dialog box that is popped up. In this mode the cursor will look like a paint brush and a dashed rectangle, the \Ibrush\i, will be attached to your cursor. As you move the brush across the plot points in the brush will be highlighted. Points outside of the brush will not be highlighted unless they are marked as selected. To select points in the brushing mode (make their highlighting permanent) hold the mouse button down as you move the brush.
\piIn the plot in Figure 10 the points within the middle of the \B\ff<(CG)Courier>hardness\b\ff<(CG)Times> range have been highlighted using a long, thin brush (you can change the size of your brush using the \BResize Brush\b command on the \BScatmat\b menu). In the plot of\B \ff<(CG)Courier>abrasion-loss\b\ff<(CG)Times> against \B\ff<(CG)Courier>tensile-strength\b\ff<(CG)Times> you can see that the highlighted points seem to follow a curve. If you want to fit a model to this data this suggests fitting a model that accounts for this curvature.\pn
\piA scatterplot matrix is also useful for examining the relationship between a quantitative variable and several categorical variables. In the data\pn
(Devore and Peck [11, page 595, Example 14]) the yield of tomato plants is recorded for an experiment run at four different planting densities and using three different varieties. In the plot in Figure 11 a long, thin brush has been used to highlight the points in the third variety. If there is no interaction between the varieties and the density then the shape of the highlighted points should move approximately in parallel as the brush is moved from one variety to another.
\piLike \B\ff<(CG)Courier>spin-plot\b\ff<(CG)Times>, the function \B\ff<(CG)Courier>scatterplot-matrix\b\ff<(CG)Times> also accepts the optional keyword argument \B\ff<(CG)Courier>scale\b\ff<(CG)Times>.\pn
\BExercises\b
\po1. Devore and Peck [11, Exercise 13.62] describe an experiment to determine the effect of oxygen concentration on fermentation end products. Four oxygen concentrations and two types of sugar were used. The data are\pn
and are on file \B\ff<(CG)Courier>oxygen.lsp\b\ff<(CG)Times>. Use a scatterplot matrix to examine these data.
2. Use scatterplot matrices to examine the data in the examples and exercises of Section 6.1 above.
\jcFigure 11: Scatterplot matrix for tomato yield data with points from
the third variety highlighted.\jf
\B6.3 Interacting with Individual Plots\b
Rotating plots and scatterplot matrices are interactive plots. Simple scatter plots also allow some interaction: If you select the \BShow Labels\b option in the plot menu a label will appear next to a highlighted point. You can use either the selecting or the brushing mode to highlight points. The default labels are of the form "0", "1", ... . (In Lisp it is conventional to start numbering indices with 0, rather than 1.) Most plotting\t<-1> \t<0>unctions allow you to supply a list of case labels using the :\B\ff<(CG)Courier>point-labels\b\ff<(CG)Times> keyword.
\piAnother option, useful in viewing large data sets, is to remove a subset of the points from your plot. This can be done by selecting the points you want to remove and then choosing \BRemove Selection\b from the plot menu. The plot can then be rescaled using the \BRescale Plot\b option. Alternatively, the \BFocus on Selection\b menu item removes all unselected points from the plot.
When a set of points is selected in a plot you can change the symbol used to display the points using the \BSelection Symbol\b item. On systems with color monitors you can set the color of selected points with the \BSelection Color\b item.
You can save the indices of the selected points in a variable by selecting the \BSelection...\b item in the plot menu. A dialog will ask you for a name for the selection. When no points are selected you can use the \BSelection...\b menu item to specify the indices of the points to select. A dialog will ask you for an expression for determining the selection to use. The expression can be any Lisp expression that evaluates to a list of indices\pn.
\B6.4 Linked Plots\b
When you brush or select in a scatterplot matrix you are looking at the interaction of a set of separate scatterplots\SP16\Sp. You can construct your own set of interacting plots by choosing the \BLink View\b option \jcFigure 12: Scatterplot and histogram with points from one sugar group highlighted.\jf
from the menus of the plots you want to link. For example, using the data from Exercise 1 in Section 6.2 we can put \B\ff<(CG)Courier>ethanol\b\ff<(CG)Times> and \B\ff<(CG)Courier>oxygen\b\ff<(CG)Times> in a scatterplot and \B\ff<(CG)Courier>sugar\b\ff<(CG)Times> in a histogram. If we link these two plots then selecting one of the two sugar groups in the histogram highlights the corresponding points in the scatterplot, as shown in Figure 12.
If you want to be able to select the points with particular labels you can use the \B\ff<(CG)Courier>name-list\b\ff<(CG)Times> function to generate a window with a list of names in it. This window can be linked with any plot, and selecting a name in a name list will then highlight the corresponding points in the linked plots. You can use the \B\ff<(CG)Courier>name-list\b\ff<(CG)Times> function with a numerical argument; e. g.
\B\ff<(CG)Courier>(name-list 10)\b\ff<(CG)Times>
will generate a list with the names "0" , ..., "9", or you can give it a list of case labels of your own.
\BExercise\b
Try to use linked scatter plots and histograms on the data in the examples and exercises of Sections 6.1 and 6.2.
\B6.5 Modifying a Scatter Plot\b
After producing a scatterplot of a data set we might like to add a line, a regression line for example, to the plot. As an example, Devore and Peck [11, page 105, Example 2] describe a data set collected to examine the effect of bicycle lanes on drivers and bicyclists. The variables given by
represent the distance between the cyclist and the roadway center line and the distance between the cyclist and a passing car, respectively, recorded in ten cases. A regression line fit to these data, with \B\ff<(CG)Courier>separation\b\ff<(CG)Times> as the dependent variable, has a slope of 0.66 and an intercept of -2.18. Let's see how to add this line to a scatterplot of the data.
to produce a scatterplot of these points. To be able to add a line to the plot, however, we must be able to refer to it within XLISP-STAT. To accomplish this let's assign the result returned by the p\B\ff<(CG)Courier>lot-points\b\ff<(CG)Times> function to a symbol\SP17\Sp:
The result returned by plot-points is an xlisp \Iobject\i. To use an object you have to send it \Imessages\i. This is done using the \B\ff<(CG)Courier>send\b\ff<(CG)Times> function, as in the expression
to tell \B\ff<(CG)Courier>myplot\b\ff<(CG)Times> to add the graph of a line\I a\i + \Ib x\i, with \Ia\i = -2.18 and \Ib\i = 0.66, to itself. The \Imessage selector\i is \B\ff<(CG)Courier>:abline\b\ff<(CG)Times>, the numbers -2.18 and 0.66 are the arguments. The message consists of the selector and the arguments. Message selectors are always Lisp keywords; that is, they are symbols that start with a colon. Before typing in this expression you might want to resize and rearrange the listener window so you can see the plot and type commands at the same time. Once you have resized the listener window you can easily switch between the small window and a full size window by clicking in the zoom box at the right corner of the window title bar. The resulting plot is shown in Figure 13.
\piScatter plot objects understand a number of other messages. One other message is the \B\ff<(CG)Courier>:help\b\ff<(CG)Times> message\SP18\Sp:\pn
\jcFigure 13: Scatterplot of bicycle data with fitted line.\jf
\piThe list of topics will be the same for all scatter plots but will be somewhat different for rotating plots, scatterplot matrices or histograms.
The \B\ff<(CG)Courier>:clear\b\ff<(CG)Times> message, as its name suggests, clears the plot and allows you to build up a new plot from scratch. Two other useful messages are \B\ff<(CG)Courier>:add-points\b\ff<(CG)Times> and \B\ff<(CG)Courier>:add-lines\b\ff<(CG)Times>. To find out how to use them we can use the \B\ff<(CG)Courier>:help\b\ff<(CG)Times> message with \B\ff<(CG)Courier>:add-points\b\ff<(CG)Times> or \B\ff<(CG)Courier>:add-lines\b\ff<(CG)Times> as arguments: \pn
\t<-1>Adds points to plot. POINTS is a list of sequences, POINT-LABELS a list of
strings. If DRAW is true the new points are added to the screen.\t<0>
> (send myplot :help :add-lines)
:ADD-LINES
Method args: (lines &key type (draw t))
\t<-1>Adds lines to plot. LINES is a list of sequences, the coordinates of the line
starts. TYPE is normal or dashed. If DRAW is true the new lines are added to the screen.\t<0>
> \b\ff<(CG)Times>\ll<-2>
\ll<1>
The plot produced above shows some curvature in the data. A regression of \B\ff<(CG)Courier>separation\b\ff<(CG)Times> on a linear and a quadratic term in \B\ff<(CG)Courier>travel-space\b\ff<(CG)Times> produces estimates of -16.41924 for the constant, 2.432667 as the coefficient of the linear term and -0.05339121 as the coefficient of the quadratic term. Let's use the \B\ff<(CG)Courier>:clear\b\ff<(CG)Times>, \B\ff<(CG)Courier>:add-points\b\ff<(CG)Times> and\B \ff<(CG)Courier>:add-lines\b\ff<(CG)Times> messages to change \B\ff<(CG)Courier>myplot\b\ff<(CG)Times> to show the data along with the fitted quadratic model. First we use the expressions
\B\ff<(CG)Courier>(def x (rseq 12 22 50))\b\ff<(CG)Times>
\jcFigure 14: Scatterplot of bicycle data with fitted curve.\jf
\B\ff<(CG)Courier>(def y (+ -16.41924 (* 2.432667 x) (* -0.05339121 (* x x))))\b\ff<(CG)Times>
to define \B\ff<(CG)Courier>x\b\ff<(CG)Times> as a grid of 50 equally spaced points between 12 and 22 and \B\ff<(CG)Courier>y\b\ff<(CG)Times> as the fitted response at these \B\ff<(CG)Courier>x\b\ff<(CG)Times> values. Then the expressions
\ll<-2>
\B\ff<(CG)Courier>(send myplot :clear)\ll<1>
(send myplot :add-points travel-space separation)
(send myplot :add-lines x y)\b\ff<(CG)Times>\ll<-2>
\ll<1>
change myplot to look like Figure 14. Of course we could have used\B \ff<(CG)Courier>plot-points\b\ff<(CG)Times> to get a new plot and just modified that plot with \B\ff<(CG)Courier>:add-lines\b\ff<(CG)Times>, but the approach used here allowed us to try out all three messages.
\B6.6 Dynamic Simulations\b
As another illustration of what you can do by modifying existing plots let's construct a dynamic simulation - a movie - to examine the variation in the shape of histograms of samples from a standard normal distribution. To start off, use the expression
\ll<-2>
\B\ff<(CG)Courier>(def h (histogram (normal-rand 20)))\b\ff<(CG)Times>
\ll<1>
to construct a single histogram and save its plot object as \B\ff<(CG)Courier>h\b\ff<(CG)Times>. The \B\ff<(CG)Courier>:clear\b\ff<(CG)Times> message is available for histograms as well. As you can see from its help message
\ll<-2>
\B\ff<(CG)Courier>> (send h :help :clear)\ll<1>
:CLEAR
Message args: (&optional (draw t))
C\t<-1>lears the plot data. If DRAW is nil the plot is redrawn; otherwise\t<0> its\b\ff<(CG)Times>
it takes an optional argument. If this argument is \B\ff<(CG)Courier>NIL\b\ff<(CG)Times> then the plot will not actually be redrawn until some other event causes it to be redrawn. This is useful for dynamic simulations. Rearrange and resize your windows until you can see the histogram window even when the listener window is the active window. Then type the expression\SP19\Sp
\ll<-2>
\B\ff<(CG)Courier>(dotimes \s(i 50)\ll<1>
\s\s(send h :clear nil)
\s\s(send h :add-points (normal-rand 20)))\b\ff<(CG)Times>\ll<-3>
\ll<1>
This should produce a sequence of 50 histograms, each based on a sample of size 20. By giving the keyword argument \B\ff<(CG)Courier>draw\b\ff<(CG)Times> with value \B\ff<(CG)Courier>NIL\b\ff<(CG)Times> to the \B\ff<(CG)Courier>:clear\b\ff<(CG)Times> message you have insured that each histogram stays on the screen until the next one is ready to be drawn. Try the example again without this argument and see what difference it makes.
\piProgrammed dynamic simulations may provide another approach to viewing the relation among several variables. As a simple example, suppose we want to examine the relation between the variables\B \ff<(CG)Courier>abrasion-loss\b\ff<(CG)Times> and \B\ff<(CG)Courier>hardness\b\ff<(CG)Times> introduced in Section 6.2 above. Let's start with a histogram of \B\ff<(CG)Courier>abrasion-loss\b\ff<(CG)Times> produced by the expression\pn
\ll<-2>
\B\ff<(CG)Courier>(def h (histogram abrasion-loss))\b\ff<(CG)Times>
\ll<1>
The messages \B\ff<(CG)Courier>:point-selected\b\ff<(CG)Times>, \B\ff<(CG)Courier>:point-hilited\b\ff<(CG)Times> and \B\ff<(CG)Courier>:point-showing\b\ff<(CG)Times> are particularly useful for dynamic simulations. Here is the help information for \B\ff<(CG)Courier>:point-selected\b\ff<(CG)Times> in a histogram:
\ll<-1>
\B\ff<(CG)Courier>> (send h :help :point-selected)\ll<1>
:POINT-SELECTED
Method args: (point &optional selected)
Sets or returns selection status (true or NIL) of POINT. Sends
:ADJUST-POINT-SCREEN-STATES message if states are set. Vectorized.
> \b\ff<(CG)Times>\ll<-2>
\ll<1>
Thus you can use this message to determine whether a point is currently selected and also to select or unselect it. Again rearrange the windows so you can see the histogram while typing in the listener window and type the expression
\s\s(send h :point-selected i nil))\b\ff<(CG)Times>\ll<-2>
\ll<1>
The expression \B\ff<(CG)Courier>(order hardness)\b\ff<(CG)Times> produces the list of indices of the ordered values of hardness. Thus the first element of the result is the index of the smallest element of hardness, the second element is the index of the second smallest element of hardness, etc.. The loop moves through each of these indices and selects and unselects the corresponding point.
\piThe result on the screen is very similar to the result of brushing a \B\ff<(CG)Courier>hardness\b\ff<(CG)Times> histogram linked to an \B\ff<(CG)Courier>abrasion-loss\b\ff<(CG)Times> histogram from left to right. The drawback of this approach is that it is harder to write an expression than to use a mouse. On the other hand, when brushing with a mouse you tend to focus your attention on the plot you are brushin\png, rather than on the other linked plots. When you run a dynamic simulation you do not have to do anything while the simulation is running and can therefore concentrate fully on the results.
An intermediate solution is possible: You can set up a dialog with a scroll bar that will run through the indices in the list \B\ff<(CG)Courier>(order hardness)\b\ff<(CG)Times>, selecting the corresponding point as it is scrolled. An example in Section 8 will show you how to do this.
Like most Lisp systems XLISP-STAT is not ideally suited to real time simulations because of the need for garbage collection, to reclaim dynamically allocated storage. This is the reason that the simulations in this section will occasionally pause for a second or two. Nevertheless, in a simple simulation like the last one each iteration may still be too fast for you to be able to pick up any pattern. To slow things down you can add some extra busy work to each iteration. For example, you could use
\B\ff<(CG)Courier>(dolist \s(i (order hardness))
\s\s(send h :point-selected i t)
\s\s(dotimes (i 100))
\s\s(send h :point-selected i nil))\b\ff<(CG)Times>
\n\ff<(CG)Times>\fs<12>\jc\ll<1>\ps<100>\t<0>\FpFigure 7: Two views of a rotatable plot of data on iron content, aluminum content and phosphate absorption in sediment samples.\jf
\jcFigure 8: Rotatable plot of measurements on permafrost samples.\jf
\SP\pi17\SpThe result returned by \B\ff<(CG)Courier>plot-points\b\ff<(CG)Times> is printed something like \B\ff<(CG)Courier>#<Object: 2010278, prototype = SCATTERPLOT-PROTO>\b\ff<(CG)Times>. This is not the value returned by the function, just its printed representation . There are several other data types that are printed this way; \Ifile streams\i, as returned by the \B\ff<(CG)Courier>open\b\ff<(CG)Times> function, are one example. For the most part you can ignore these printed results. There is one unfortunate feature, however: the form \B\ff<(CG)Courier>#<...>\b\ff<(CG)Times> means that there is no printed form of this data type that the Lisp reader can understand. As a result, if you forget to give your plot a name you can't cut and paste the result into a \B\ff<(CG)Courier>def\b\ff<(CG)Times> expression -- you have to redo the plot or use the history mechanism.
\SP18\Sp\pnTo keep things simple I will use the term \Imessage \ito refer to a message corresponding to a message selector.
\SP\ff<(CG)Times>\pi19\Sp\ff<(CG)Courier>dotimes\b\ff<(CG)Times> is one of several Lisp looping constructs. It isa special form with the syntax \B\ff<(CG)Courier>(dotimes (var count) expr ....)\b\ff<(CG)Times>. The loop is repeate\pnd \B\ff<(CG)Courier>count\b\ff<(CG)Times> times, with \B\ff<(CG)Courier>var\b\ff<(CG)Times> bound to 0, 1, ..., \B\ff<(CG)Courier>count\b\ff<(CG)Times> - 1. Other looping constructs are \B\ff<(CG)Courier>dolist\b\ff<(CG)Times>, \B\ff<(CG)Courier>do\b\ff<(CG)Times> and \B\ff<(CG)Courier>do*\b\ff<(CG)Times>.
\SP\pi20\SpRecall from Section 6.5 that \B\ff<(CG)Courier>#<Object: 1966006, prototype = REGRESSION-MODEL-PROTO>\b\ff<(CG)Times> is the printed representation of the model object returned by \B\ff<(CG)Courier>regression-model\b\ff<(CG)Times>. Unfortunately you can't cut and paste it\pn into the \B\ff<(CG)Courier>def\b\ff<(CG)Times>, but of course you can cut and paste the \B\ff<(CG)Courier>regression-model\b\ff<(CG)Times> expression or use the history mechanism.
Regression models have been implemented using XLISP-STAT's object and message sending facilities. These were introduced above in Section 6.5. You might want to review that section briefly before reading on.
\piLet's fit a simple regression model to the bicycle data of Section 6.5. The dependent variable is \B\ff<(CG)Courier>separation\b\ff<(CG)Times> and the independent variable is \B\ff<(CG)Courier>travel-space\b\ff<(CG)Times>. To form a regression model use the \B\ff<(CG)Courier>regression-model\b\ff<(CG)Times> function:\pn
The basic syntax for the \B\ff<(CG)Courier>regression-model\b\ff<(CG)Times> function is
\B\ff<(CG)Courier>(regression-model x y)\b\ff<(CG)Times>
For a simple regression \B\ff<(CG)Courier>x\b\ff<(CG)Times> can be a single list or vector. For a multiple regression \B\ff<(CG)Courier>x\b\ff<(CG)Times> can be a list of lists or vectors or a matrix. The \B\ff<(CG)Courier>regression-model\b\ff<(CG)Times> function also takes several optional keyword arguments. The most important ones are \B\ff<(CG)Courier>:intercept\b\ff<(CG)Times>, \B\ff<(CG)Courier>:print\b\ff<(CG)Times>, and \B\ff<(CG)Courier>:weights\b\ff<(CG)Times>. Both \B\ff<(CG)Courier>:intercept\b\ff<(CG)Times> and \B\ff<(CG)Courier>:print\b\ff<(CG)Times> are \B\ff<(CG)Courier>T\b\ff<(CG)Times> by default. To get a model without an intercept use the expression
\ll<-2>
\B\ff<(CG)Courier>(regression-model x y :intercept nil)\b\ff<(CG)Times>
\ll<1>
To form a weighted regression model use the expression
\ll<-2>
\B\ff<(CG)Courier>(regression-model x y :weights w)\b\ff<(CG)Times>
\ll<1>
where \B\ff<(CG)Courier>w\b\ff<(CG)Times> is a list or vector of weights the same length as \B\ff<(CG)Courier>y\b\ff<(CG)Times>. The variances of the errors are assumed to be inversely proportional to the weights \B\ff<(CG)Courier>w\b\ff<(CG)Times>.
\piThe \B\ff<(CG)Courier>regression-model\b\ff<(CG)Times> function prints a very simple summary of the fit model and returns a model object as its result. To be able to examine the model further assign the returned model object to a variable using an expression like\SP20\Sp\pn
I have given the keyword argument \B\ff<(CG)Courier>:print nil\b\ff<(CG)Times> to suppress the printing of the summary, since we have already seen it. To find out what messages are available use the \B\ff<(CG)Courier>:help\b\ff<(CG)Times> message:
Many of these messages are self explanatory, and many have already been used by the \B\ff<(CG)Courier>:display\b\ff<(CG)Times> message, which \B\ff<(CG)Courier>regression-model\b\ff<(CG)Times> sends to the new model to print the summary. As examples let's try the \B\ff<(CG)Courier>:coef-estimates\b\ff<(CG)Times> and \B\ff<(CG)Courier>:coef-standard-errors\b\ff<(CG)Times> messages\SP2\Sp:
The \B\ff<(CG)Courier>:plot-residuals\b\ff<(CG)Times> message will produce a residual plot. To find out what residuals are plotted against let's look at the help information:
we can construct two plots of the data as shown in Figure 15. By linking the plots we can use the mouse to identify points in both plots simultaneously. A point that stands out is observation 6 (starting the count at 0, as usual).
\jcFigure 15: Linked raw data and residual plots for the bicycles example.\jf
The plots both suggest that there is some curvature in the data; this curvature is particularly pronounced in the residual plot if you ignore observation 6 for the moment. To allow for this curvature we might try to fit a model with a quadratic term in \B\ff<(CG)Courier>travel-space\b\ff<(CG)Times>:
I have used the exponentiation function "^" to compute the square of \B\ff<(CG)Courier>travel-space\b\ff<(CG)Times>. Since I am now forming a multiple regression model the first argument to \B\ff<(CG)Courier>regression-model\b\ff<(CG)Times> is a list of the \B\ff<(CG)Courier>x\b\ff<(CG)Times> variables.
\piYou can proceed in many directions from this point. If you want to calculate Cook's distances for the observations you can first compute internally studentized residuals as\pn
The seventh entry - observation 6, counting from zero - clearly stands out.
\piAnother approach to examining residuals for possible outliers is to use the Bayesian residual plot proposed by Chaloner and Brant [7], which can be obtained using the message :\B\ff<(CG)Courier>plot-bayes-residuals\b\ff<(CG)Times>. The expression \B\ff<(CG)Courier>(send bikes2 :plot-bayes-residuals)\b\ff<(CG)Times> produces the plot in Figure 16.\pn The bars represent mean
\ff<(CG)Symbol>
\ff<(CG)Times>2SD of the posterior distribution of the actual realized errors, based on an improper uniform prior distribution on the regression coefficients. The \Iy\i axis is in units of \ff<(CG)Symbol>s\ff<(CG)Times>. Thus this plot suggests the probability that point 6 is three or more standard deviations from the mean is about 3%; the probability that it is at least two standard deviations from the mean is around 50%.
\BExercises\b
\po1. Using the variables \B\ff<(CG)Courier>absorption\b\ff<(CG)Times>, \B\ff<(CG)Courier>iron\b\ff<(CG)Times> and \B\ff<(CG)Courier>aluminum\b\ff<(CG)Times> introduced in Section 6.1 above construct and examine a model with \B\ff<(CG)Courier>absorption\b\ff<(CG)Times> as the dependent variable.
2. Using the variables \B\ff<(CG)Courier>abrasion-loss\b\ff<(CG)Times>, \B\ff<(CG)Courier>tensile-strength\b\ff<(CG)Times> and \B\ff<(CG)Courier>hardness\b\ff<(CG)Times> introduced in Section 6.2 above construct and examine a model with \B\ff<(CG)Courier>abrasion-loss\b\ff<(CG)Times> as the dependent variable.\pn
\SP\pi21\Sp\fs<10>Ordinarily the entries in the lists returned by these messages correspond simply to the intercept, if one is included in the model, followed by the independent variables as they were supplied to \B\ff<(CG)Courier>regression-model\b\ff<(CG)Times>. However, if degeneracy is detected during computations some variables will not be used in the fit; they will be marked as \Ialiased\i in the printed summary. The indices of the variables used can be obtained by the :\ff<(CG)Courier>basis \ff<(CG)Times>message; the entries in the list returned by :\B\ff<(CG)Courier>coef-estimates\b\ff<(CG)Times> correspond to the intercept, if appropriate, followed by the coefficients of the elements in the basis. The messages \B\ff<(CG)Courier>:x-matrix\b\ff<(CG)Times> and \B\ff<(CG)Courier>:xtxinv\b\ff<(CG)Times> are similar in that they use only the variables in the basis.
\SP\pi22\SpThe / function is used here with three arguments. The first argument is divided by the second, and the result is then divided by the third. Thus the result of the expression (/ 6 3 2) is 1.\pn
C\n\B\ff<(CG)Times>\fs<14>\jf\ll<1>\ps<100>\t<0>\Fp8 Defining Your Own Functions and Methods\b\fs<12>
This section gives a brief introduction on how to program XLISP-STAT. The most basic programming operation is to define a new function. Closely related is the idea of defining a new method for an object.\SP23\Sp
\B8.1 Defining Functions\b
You can use the XLISP language to define functions of your own. Many of the functions you have been using so far are written in this language. The special form used for defining functions is called \B\ff<(CG)Courier>defun\b\ff<(CG)Times>. The simplest form of the \B\ff<(CG)Courier>defun\b\ff<(CG)Times> syntax is
\B\ff<(CG)Courier>(defun fun args expression)\b\ff<(CG)Times>
where \B\ff<(CG)Courier>fun\b\ff<(CG)Times> is the symbol you want to use as the function name, \B\ff<(CG)Courier>args\b\ff<(CG)Times> is the list of the symbols you want to use as arguments and \B\ff<(CG)Courier>expression\b\ff<(CG)Times> is the body of the function. Suppose for example that you want to define a function to delete a case from a list. This function should take as its arguments the list and the index of the case you want to delete. The body of the function can be based on either of the two approaches described in Section 4.3 above. Here is one approach:
\B\ff<(CG)Courier>(defun delete-case (x i)
\s(select x (remove i (iseq 0 (- (length x) 1))) ) )\b\ff<(CG)Times>
I have used the function \B\ff<(CG)Courier>length\b\ff<(CG)Times> in this definition to determine the length of the argument \B\ff<(CG)Courier>x\b\ff<(CG)Times>. Note that none of the arguments to \B\ff<(CG)Courier>defun\b\ff<(CG)Times> are quoted: \B\ff<(CG)Courier>defun\b\ff<(CG)Times> is a special form that does not evaluate its arguments.
\piUnless the functions you define are very simple you will probably want to define them in a file and load the file into XLISP-STAT with the \B\ff<(CG)Courier>load\b\ff<(CG)Times> command. You can put the functions in the \B\ff<(CG)Courier>statinit.lsp\b\ff<(CG)Times> file or include in \B\ff<(CG)Courier>statinit.lsp \b\ff<(CG)Times>a \B\ff<(CG)Courier>load\b\ff<(CG)Times> command that will load another file. The version of XLISP-STAT for the Macintosh includes a simple editor that can be used from within XLISP-STAT. The editor is described briefly in Section 5.6 above.
You can also write functions that send messages to objects. Here is a function that takes two regression models, assumed to be nested, and \pncomputes the \IF\i statistic for comparing these models:
\B\ff<(CG)Courier>(defun f-statistic (m1 m2)
"Args: (m1 m2)
Computes the F statistic for testing model m1 within model m2."
This definition uses the Lisp \B\ff<(CG)Courier>let\b\ff<(CG)Times> construct to establish some local \Ivariable bindings\i. The variables \B\ff<(CG)Courier>ss1\b\ff<(CG)Times>, \B\ff<(CG)Courier>df1\b\ff<(CG)Times>, \B\ff<(CG)Courier>ss2\b\ff<(CG)Times> and \B\ff<(CG)Courier>df2\b\ff<(CG)Times> are defined in terms of the two model arguments \B\ff<(CG)Courier>m1\b\ff<(CG)Times> and \B\ff<(CG)Courier>m2\b\ff<(CG)Times>, and are then used to compute the \IF\i statistic. The string following the argument list is a \Idocumentation string\i. When a documentation string is present in a \B\ff<(CG)Courier>defun\b\ff<(CG)Times> expression \B\ff<(CG)Courier>defun\b\ff<(CG)Times> will install it so the \B\ff<(CG)Courier>help\b\ff<(CG)Times> function will be able to retrieve it.
\B8.2 Anonymous Functions\b
Suppose you would like to plot the function \If\i(\Ix\i) = 2\Ix\i + \Ix\i\SP2\Sp over the range -2 \ff<(CG)Symbol>
\ff<(CG)Times> \Ix\i \ff<(CG)Symbol>
\ff<(CG)Times> 3. We can do this by first defining a function \B\ff<(CG)Courier>f\b\ff<(CG)Times> and then using \B\ff<(CG)Courier>plot-function\b\ff<(CG)Times>:
\ll<-2>
\B\ff<(CG)Courier>(defun f (x) (+ (* 2 x) (^ x 2)))\ll<1>
(plot-function #'f -2 3)\b\ff<(CG)Times>\ll<-2>
\ll<1>
This is not too hard, but it nevertheless involves one unnecessary step: naming the function \B\ff<(CG)Courier>f\b\ff<(CG)Times>. You probably won't need this function again; it is a throw away function defined only to allow you to give it to \B\ff<(CG)Courier>plot-function\b\ff<(CG)Times> as an argument. It would be nice if you could just give \B\ff<(CG)Courier>plot-function\b\ff<(CG)Times> an expression that constructs the function you want. Here is such an expression:
\ll<-2>
\B\ff<(CG)Courier>(function (lambda (x) (+ (* 2 x) (^ x 2))))\b\ff<(CG)Times>
\ll<1>
There are two steps involved. The first is to specify the definition of your function. This is done using a\I lambda expression\i, in this case
\ll<-2>
\B\ff<(CG)Courier>(lambda (x) (+ (* 2 x) (^ x 2)))\b\ff<(CG)Times>
\ll<1>
A lambda expression is a list starting with the symbol \B\ff<(CG)Courier>lambda\b\ff<(CG)Times>, followed by the list of arguments and the expressions making up the body of the function. The lambda expression is only a definition, it is not yet a function, a piece of code that can be applied to arguments. The special form \B\ff<(CG)Courier>function\b\ff<(CG)Times> takes the lambda list and constructs such a function. The result can be saved in a variable or it can be passed on to another function as an argument. For our plotting problem you can use the single expression
Since the function used in these expressions is never named it is sometimes called an \Ianonymous function\i.
You can also construct a rotating plot of a function of two variables using the function \B\ff<(CG)Courier>spin-function\b\ff<(CG)Times>. As an example, the expression
\ll<-2>
\B\ff<(CG)Courier>(spin-function #'(lambda (x y) (+ (^ x 2) (^ y 2))) -1 1 -1 1)\b\ff<(CG)Times>
\ll<1>
constructs a plot of the function \If\i(\Ix\i,\I y\i) = \Ix\i\SP2\Sp + \Iy\i\SP2\Sp over the range -1\ff<(CG)Symbol>
\ff<(CG)Times> \Ix\i \ff<(CG)Symbol>
\ff<(CG)Times> 1, -1\ff<(CG)Symbol>
\ff<(CG)Times> \Iy\i \ff<(CG)Symbol>
\ff<(CG)Times> 1 using a 6 x 6 grid.
\piThe number of grid points can be changed using the \B\ff<(CG)Courier>:num-points\b\ff<(CG)Times> keyword. The result is shown in Figure 17. Again it convenient to use a lambda expression to specify the function to be plotted. There are a number of situations in which you might want to pass a function on as an argument without first going through the trouble of thinking up a name for the function and defining it using defun. A few additional examples are given in the next subsection.\pn
\B8.3 Some Dynamic Simulations\b
In Section 6.6 we used a loop to control a dynamic simulation in which points in a histogram of one variable were selected and deselected in the order of a second variable. Let's look at how to run the same simulation using a \Islider\i to control the simulation.
\piA slider is a modeless dialog box containing a scroll bar and a value display field. As the scroll bar is\pn moved the displayed value is changed and an action is taken. The action is defined by an \Iaction function\i given to the scroll bar, a function that is called with one value, the current slider value, each time the value is changed by the user. There are two kinds of sliders, sequence sliders and
interval sliders. Sequence sliders take a sequence (a list or a vector) and scroll up and down the sequence. The displayed value is either the current sequence element or the corresponding element of a display sequence. An interval slider dialog takes an interval, divides it into a grid and scrolls through the grid. By default a grid of around 30 points is used; the exact number and the interval end points are adjusted to give nice printed values. The current interval point is displayed in the display field.
\piFor our example lets use a sequence slider to scroll through the elements of the \B\ff<(CG)Courier>hardness\b\ff<(CG)Times> list in order and select the corresponding element of \B\ff<(CG)Courier>abrasion-loss\b\ff<(CG)Times>. The expression \pn
\ll<-2>
\B\ff<(CG)Courier>(def h (histogram abrasion-loss))\b\ff<(CG)Times>
\ll<1>
sets up a histogram and saves its plot object in the variable \B\ff<(CG)Courier>h\b\ff<(CG)Times>. The function \B\ff<(CG)Courier>sequence-slider-dialog\b\ff<(CG)Times> takes a list or vector and opens a sequence slider to scroll through its argument. To do something useful with this dialog we need to give it an action function as the value of the keyword argument \B\ff<(CG)Courier>:action\b\ff<(CG)Times>. The function should take one argument, the current value of the sequence controlled by the slider. The expression
\s\s\s\s\s\s(send h :point-selected i t)))\b\ff<(CG)Times>\ll<-2>
\ll<1>
sets up a slider for moving the selected point in the \B\ff<(CG)Courier>abrasion-loss\ff<(CG)Times> \bhistogram along the order of the \B\ff<(CG)Courier>hardness\b\ff<(CG)Times> variable. The histogram and scrollbar are shown in Figure 18. The action function is called every time the slider is moved. It is called with the current element of the sequence \B\ff<(CG)Courier>(order hardness)\b\ff<(CG)Times>, the index of the point to select. It clears all current selections and then selects the point specified in the call from the slider. The body of the function is almost identical to the body of the \B\ff<(CG)Courier>dotimes\b\ff<(CG)Times> loop used in Section 6.6. The slider is thus an interactive, graphically controlled version of this loop.
\jcFigure 18: Slider-controlled animation of a histogram.\jf
\piAs another example, suppose we would like to examine the effect of changing the exponent in a Box-Cox power transformation\pn
on a probability plot. As a first step we might define a function to compute the power transformation and normalize the result to fall between zero and one:
This definition uses the \B\ff<(CG)Courier>let*\b\ff<(CG)Times> form to establish some local variable bindings. The \B\ff<(CG)Courier>let*\b\ff<(CG)Times> form is like the \B\ff<(CG)Courier>let\b\ff<(CG)Times> form used above except that \B\ff<(CG)Courier>let*\b\ff<(CG)Times> defines its variables sequentially, allowing a variable to be defined in terms of other variables already defined in the \B\ff<(CG)Courier>let*\b\ff<(CG)Times> expression. In this case the variables \B\ff<(CG)Courier>min\b\ff<(CG)Times> and \B\ff<(CG)Courier>max\b\ff<(CG)Times> are defined in terms of the variable \B\ff<(CG)Courier>bcx\b\ff<(CG)Times>.
\piNext we need a set of positive numbers to transform. Let's use a sample of thirty observations from a \ff<(CG)Symbol>c\SB\t<-1>4\Sb\ff<(CG)Times>\t<-4>2\t<0> distribution and order the observations\pn
\B\ff<(CG)Courier>(def x (sort-data (chisq-rand 30 4)))\b\ff<(CG)Times>
The normal quantiles of the expected uniform order statistics are given by
\B\ff<(CG)Courier>(def r (normal-quant (/ (iseq 1 30) 31)))\b\ff<(CG)Times>
and a probability plot of the untransformed data is constructed using
\ll<-2>
\B\ff<(CG)Courier>(def myplot (plot-points r (bc x 1)))\b\ff<(CG)Times>
\ll<1>
Since the power used is 1 the function \B\ff<(CG)Courier>bc\b\ff<(CG)Times> just rescales the data.
\piThere are several ways to set up a slider dialog to control the power parameter. The simplest approach is to use the function \B\ff<(CG)Courier>interval-slider-dialog\b\ff<(CG)Times>:\pn
\s\s\s\s\s\s(send\t<-1> \t<0>myplot\t<-2> \t<0>:add-points r (bc x p))))\b\ff<(CG)Times>
\B\ff<(CG)Courier>interval-slider-dialog\b\ff<(CG)Times> takes a list of two numbers, the lower and upper bound of an interval. The action function is called with the current position in the interval.
\piThis approach works fine on a Mac II but may be a bit slow on a Mac Plus or a Mac SE. An alternative is to precompute the transformations for a list of powers and then use the precomputed values in the display. For example, using the powers defined in\pn
we can compute the transformed data for each power and save the result as the variable \B\ff<(CG)Courier>xlist\b\ff<(CG)Times>:
\B\ff<(CG)Courier>(def xlist (mapcar #'(lambda (p) (bc x p)) powers))\b\ff<(CG)Times>
The function \B\ff<(CG)Courier>mapcar\b\ff<(CG)Times> is one of several \Imapping functions \iavailable in Lisp. The first argument to \B\ff<(CG)Courier>mapcar\b\ff<(CG)Times> is a function. The second argument is a list. \B\ff<(CG)Courier>mapcar\b\ff<(CG)Times> takes the function, applies it to each element of the list and returns the list of the results\SP24\Sp. Now we can use a sequence slider to move up and down the transformed values in \B\ff<(CG)Courier>xlist\b\ff<(CG)Times>:
\s\s\s\s\s\s(send myplot :add-points r x)))\b\ff<(CG)Times>\ll<-2>
\ll<1>
Note that we are scrolling through a list of lists and the element passed to the action function is thus the list of current transformed values. We would not want to see these values in the display field on the slider, so I have used the keyword argument \B\ff<(CG)Courier>:display\b\ff<(CG)Times> to specify an alternative display sequence, the powers used in the transformation.
\B8.4 Defining Methods\b
When a message is sent to an object the object system will use the object and the method selector to find the appropriate piece of code to execute. Different objects may thus respond differently to the same message. A linear regression model and a nonlinear regression model might both respond to a \B\ff<(CG)Courier>:coef-estimates\b\ff<(CG)Times> message but they will execute different code to compute their responses.
\piThe code used by an object to respond to a message is called a \Imethod\i. Objects are organized in a hierarchy in which objects inherit from other objects. If an object does not have a met\pnhod of its own for responding to a message it will use a method inherited from one of its ancestors. The \B\ff<(CG)Courier>send\b\ff<(CG)Times> function will move up the \Ipredecence list \iof an object's ancestors until a method for a message is found.
\piMost of the objects encountered so far inherit directly from \Iprototype\i objects. Scatterplots inherit from \B\ff<(CG)Courier>scatterplot-proto\b\ff<(CG)Times>, histograms from \B\ff<(CG)Courier>histogram-proto\b\ff<(CG)Times>, regression models from \B\ff<(CG)Courier>regression-model-proto\b\ff<(CG)Times>. Prototypes are just like any other object. They are essentially \Itypical\i versions of a certain kind of object that define default behavior. Almost all methods are owned by prototypes. But any object can own a method, and in the process of debugging a new method it is often better to attach the method to a separate object constructed for that purpose instead of the prototype.
To add a method to an object you can use the \B\ff<(CG)Courier>defmeth\b\ff<(CG)Times> macro. As an example, in Section 7 we calculated Cook's distances for a regression model. If you find that you are doing this very frequently then you might want to define a \B\ff<(CG)Courier>:cooks-distances\b\ff<(CG)Times> method. The general form of a method definition is:\pn
\B\ff<(CG)Courier>object\b\ff<(CG)Times> is the object that is to own the new method. In the case of regression models this can be either a specific regression model or the prototype that all regression models inherit from, \B\ff<(CG)Courier>regression-model-proto\b\ff<(CG)Times>. The argument \B\ff<(CG)Courier>:new-method\b\ff<(CG)Times> is the message selector for your method; in our case this would be \B\ff<(CG)Courier>:cooks-distances\b\ff<(CG)Times>. The argument \B\ff<(CG)Courier>arg-list\b\ff<(CG)Times> is the list of argument names for your method, and \B\ff<(CG)Courier>body\b\ff<(CG)Times> is one or more expressions making up the body of the method. When the method is used each of these expressions will be evaluated, in the order in which they appear.
\piHere is an expression that will install the\B \ff<(CG)Courier>:cooks-distances\b\ff<(CG)Times> method:\pn
The variable \B\ff<(CG)Courier>self\b\ff<(CG)Times> refers to the object that is receiving the message.
\B8.5 Plot Methods\b
All plot activities are accomplished by sending messages to plot objects. For example, every time a plot needs to be redrawn the system sends the plot the \B\ff<(CG)Courier>:redraw\b\ff<(CG)Times> message. By defining a new method for this message you can change the way a plot is drawn. Similarly, when the mouse is moved or clicked in a plot the plot is sent the \B\ff<(CG)Courier>:do-mouse\b\ff<(CG)Times> message. Menu items also send messages to their plots when they are selected. If you are interested in modifying plot behavior you may be able to get started by looking at the methods defined in the graphics files loaded on start up. Further details will be given in [17].
\B\fs<14>9 Matrices and Arrays\b\fs<12>
XLISP-STAT includes support for multidimensional arrays patterned after the Common Lisp standard described in detail in Steele [15]. The funcitons supported are listed in Section C.6 of the appendix.
\piIn addition to the standard Common Lisp array functions XLISP-STAT also includes a number of linear algebra functions such as \B\ff<(CG)Courier>inverse\b\ff<(CG)Times>, \B\ff<(CG)Courier>solve\b\ff<(CG)Times>, \B\ff<(CG)Courier>transpose\b\ff<(CG)Times>, c\B\ff<(CG)Courier>hol-decomp\ff<(CG)Times>,\b l\B\ff<(CG)Courier>u-decomp\b\ff<(CG)Times> and \B\ff<(CG)Courier>sv-decomp\b\ff<(CG)Times>. These functions are list in the appendix as well.
An array is printed using the standard Common Lisp format. For example, a 2 by 3 matrix with rows (1 2 3) and (4 5 6) is printed as
The prefix \B\ff<(CG)Courier>#2A\b\ff<(CG)Times> indicates that this is a two-dimensional array. This form is not particularly readable, but it has the advantage that it can be pasted into expressions and will be read as an array by the XLISP reader.\SP25\Sp For matrices you can use the function \B\ff<(CG)Courier>print-matrix\b\ff<(CG)Times> to get a slightly more readable representation:
\piThe \B\ff<(CG)Courier>select\b\ff<(CG)Times> funciton can be used to extract elements or subarrays from an array. If \B\ff<(CG)Courier>A\b\ff<(CG)Times> is a two dimensional array then the expression
\B\ff<(CG)Courier>\pn(select a 0 1)\b\ff<(CG)Times>
will return element 1 of row 0 of \B\ff<(CG)Courier>A\b\ff<(CG)Times>. The expression
\B\ff<(CG)Courier>(select a (list 0 1)\M(list 0 1))\b\ff<(CG)Times>
returns the upper left hand corner of \B\ff<(CG)Courier>A\b\ff<(CG)Times>.\pi
\B\fs<14>10 Nonlinear Regression\b\fs<12>
XLISP-STAT allows the construction of nonlinear, normal regression models. The present implementation is experimental (even more experimental than the rest of the program). The definitions needed for nonlinear regression are in the file \B\ff<(CG)Courier>nonlin.lsp\b\ff<(CG)Times> on the distribution disk. This file is not loaded automatically at start up; you should load it now, using the \BLoad\b item on the \BFile\b menu or the \B\ff<(CG)Courier>load\b\ff<(CG)Times> command, to carry out the calculations in this section.
\piAs an example,Bates and Watts [1, A1.3] descrone an experiment on the relation between the velocity of an enzymatic reaction, \Iy\i, and the substrate concentration, \Ix\i. The data for an experiment in which the enzyme was treated with Puromycin are given by\pn
often provides a good model for the dependence of velocity on substrate concentration. Assuming the Michaelis-Menton function as the mean velocity at a given concentration level the function \B\ff<(CG)Courier>f1\b\ff<(CG)Times> defined by
\ll<-2>
\B\ff<(CG)Courier>(defun f1 (b) (/ (* (select b 0) x) (+ (select b 1) x)))\b\ff<(CG)Times>
\ll<1>
computes the list of mean response values at the points in \B\ff<(CG)Courier>x1\ff<(CG)Times> \bfor a parameter list \B\ff<(CG)Courier>b\b\ff<(CG)Times>. Using these definitions, which are contained in the file \B\ff<(CG)Courier>puromycin.lsp\b\ff<(CG)Times> in the \B\ff<(CG)Courier>Data\b\ff<(CG)Times> folder of the distribution disk, we can construct a nonlinear regression model using the \B\ff<(CG)Courier>nreg-model\b\ff<(CG)Times> function.
\piFirst we need initial estimates for the two model parameters. Examining the expression for the Michaelis-Menton mdoel shows that as \Ix\i increases the function approaches an asymptote, \ff<(CG)Symbol>q\SB\ff<(CG)Times>1\Sb. The second parameter, \ff<(CG)Symbol>q\SB\ff<(CG)Times>2\Sb, can be interpolated as the value of \Ix\i at which the function has reached half its asymptotic value. Using these interpretations for the parameters and a plot constructed by the expression
\pn\ll<-2>
\B\ff<(CG)Courier>(plot-points x1 y1)
\ll<1>
\b\ff<(CG)Times>shown in Figure 19 we can read off reasonable initial estimates of 200 for \ff<(CG)Symbol>q\SB\ff<(CG)Times>1\Sb and 0.1 for \ff<(CG)Symbol>q\SB\ff<(CG)Times>2\Sb. The nreg-model function takes the mean function, the response vector and an initial estimate of the parameters, computes more accurate estimates using an iteractive algorithm starting from the initial guess, and prints a summary of the results. It returns a nonlinear regression model object:\SP26\Sp
\b\ff<(CG)Times>\jcFigure 19: Plot of reaction veloctiy against substrate concentration for Puromycin experiment
\B\ff<(CG)Courier>\jf
Parameter 0: 212.684 (6.94715)
Parameter 1: 0.0641213 (0.00828094)
R Squared: 0.961261
Sigma hat: 0.004483935
Number of cases: 12
Degrees of freedom: 10
PUROMYCIN
>\b\ff<(CG)Times>
The function \B\ff<(CG)Courier>nreg-model\b\ff<(CG)Times> also takes several keyword arguments, including \B\ff<(CG)Courier>:epsilon\b\ff<(CG)Times> to specify a convergence criterion and \B\ff<(CG)Courier>:count-limit\b\ff<(CG)Times>, a limit on the number of iterations. By default these values are .001 and 20, respectively. The algorithm for fitting the model is a simple Gauss-Newton algorithm with backtracking; derivatives are computed numerically.
\piTo see how you can analyze the model further you can send \B\ff<(CG)Courier>puromycin\b\ff<(CG)Times> the \B\ff<(CG)Courier>:help\b\ff<(CG)Times> message. The result is very similar to the help information for a linear regression model. The reason for this is simple: nonlinear regression models have been implemented as objects, with the nonlinear regression model prototype \B\ff<(CG)Courier>nreg-model-proto\b\ff<(CG)Times> inheriting from of the linear regression model prototype. The internal data, the method for computing estimates and the method of computing fitted values have been modified, but the other methods remain unchanged. Once the model has been fit the Jacobian of the mean function at the estimated parameter values is used as the \IX\i matrix. In terms of this definition most of the methods for linear regression, such as the methods \B\ff<(CG)Courier>:coef-standard-errors\b\ff<(CG)Times> and \B\ff<(CG)Courier>:leverages\b\ff<(CG)Times>, still make sense, at least as first order asymptotic approximations.\pn
\piIn addition to the messages for linear regression models a nonlinear regression model can respond to the messages
\B\ff<(CG)Courier>\pn:COUNT-LIMIT
:EPSILON
:MEAN-FUNCTION
:NEW-INITIAL-GUESS
:PARAMETER-NAMES
:THETA-HAT
:VERBOSE\b\ff<(CG)Times>
\BExercices
1. Examine the residuals of the \B\ff<(CG)Courier>puromycin\b\ff<(CG)Times> model.
\po2. The file \B\ff<(CG)Courier>puromycin.lsp\b\ff<(CG)Times> also contains data \B\ff<(CG)Courier>x2\b\ff<(CG)Times> and \B\ff<(CG)Courier>y2\b\ff<(CG)Times> and mean function \B\ff<(CG)Courier>f2\b\ff<(CG)Times> for an experiment run without the Puromycin treatment. Fit a model to this data and compare the results to the experiment with Puromycin.\pn
\B\fs<14>11 One Way ANOVA\b\fs<12>
XLISP-STAT allows the construction of normal one way analysis of variance models. The definitions needed are in the file \B\ff<(CG)Courier>oneway.lsp\b\ff<(CG)Times> on the distribution disk. This file is not loaded automatically at start up; you should load it now, using the \BLoad\b item on the \BFile\b menu or the \B\ff<(CG)Courier>load\b\ff<(CG)Times> command, to carry out the calculations in this section.
As an example, suppose we would like to model the data on cholesterol levels in rural and urban Guatemalans, examined in Section 3.2, as a one way ANOVA model. The boxplots we obtained in Section 3.2 showed that the samples were skewed and the center and spread of the urban sample were larger than the center and spread of the rural sample. To compensate for these facts I will use a normal ANOVA model for the logarithms of the data:
The function \B\ff<(CG)Courier>oneway-model\b\ff<(CG)Times> requires one argument, a list of the lists or vectors representing the samples for the different groups. The model \B\ff<(CG)Courier>cholesterol\b\ff<(CG)Times> can respond to all regression messages as well as a few new ones. The new ones are
\B\ff<(CG)Courier>:BOXPLOTS
:ERROR-MEAN-SQUARE
:ERROR-DF
:GROUP-DF
:GROUP-MEAN-SQUARE
:GROUP-NAMES
:GROUP-SUM-OF-SQUARES
:GROUPED-DATA
:STANDARD-DEVIATIONS\b\ff<(CG)Times>
\B\fs<14>12 Maximization and Maximum Likelihood Estimation \b\fs<12>
XLISP-STAT includes two functions for maximizing functions of several variables. The definitions needed are in the file \B\ff<(CG)Courier>maximize.lsp\b\ff<(CG)Times> on the distribution disk. This file is not loaded automatically at start up; you should load it now, using the \BLoad\b item on the \BFile\b menu or the \B\ff<(CG)Courier>load\b\ff<(CG)Times> command, to carry out the calculations in this section. The material in this seciton is somewhat more advanced as it assumes you are familiar with the basic concepts of maximum likelihood estimation.
\piTwo functions are available for maximization. The first is \B\ff<(CG)Courier>newtonmax\b\ff<(CG)Times>, which takes a function and a list or vector representing an initial guess for the location of the maximum and attempts to find the maximum using an algorithm based on Newton's method with backtracking. The algorithm is based on the unconstrained minimization system described in Dennis and Schnabel [10].
As an example, Cox and Snell [9, Example T] describe data collected on times (in operating hours) between failures of air conditioning units on several aircraft. The data for one of the aircraft can be enterdd as
A simple model for these data might be to assume the times between failures are independent random variables with a common gamma distribution. The density of the gamma distribution can be written as
where \ff<(CG)Symbol>m\ff<(CG)Times> is the mean time between failures and \ff<(CG)Symbol>\ll<-7>b\ff<(CG)Times>\ll<1> is the gamma shape parameter. The log likelihood for the sample is thus given by\ll<-1>
We can define a Lisp function to evaluate this log likelihood by
\ll<-2>
\B\ff<(CG)Courier>(defun gllik (theta)\ll<1>
\T (let* ((mu (select theta 0))
\s(beta (select theta 1))
\s(n (length x))
\s(byn (* x (/ beta mu))))
\T\T\T\T(+ (* n (- (log beata)\M(log mu)\M(log-gamma beta)))
\s(sum (*\M(- beta 1)\M(log bym)))
\s(sum (- bym)))))\b\ff<(CG)Times>\ll<-2>
\ll<1>
This definition uses the function \B\ff<(CG)Courier>log-gamma\b\ff<(CG)Times> to evaluate \ff<(CG)Symbol>\ll<-7>G\ff<(CG)Times>\ll<1>(\ff<(CG)Symbol>\ll<-7>b\ff<(CG)Times>\ll<1>). The data and function definitions are contained in the file \B\ff<(CG)Courier>aircraft.lsp\b\ff<(CG)Times> in the \B\ff<(CG)Courier>Data\b\ff<(CG)Times> folder of the distribution disk.
\piClosed form maximum likelhood estimates are not available for the shape parameter of this distribution, but we can use \B\ff<(CG)Courier>newtonmax\b\ff<(CG)Times> to compute estimates numerically.\SP27\Sp We need an initial guess to use as a starting value in the maximization. To get initial estimates we can compute the mean and standard deviation of \B\ff<(CG)Courier>x\b\ff<(CG)Times>\pn
\ll<-2>
\B\ff<(CG)Courier>> (mean x)\ll<1>
83.5172
> (standard-deviation x)
70.8059\b\ff<(CG)Times>
\s\s\s\s\s \SB\ll<-5>^\s\s ^ ^\Sb
and use the method of moments estimates \t<-3> \ff<(CG)Symbol>\t<-1>m\ff<(CG)Times>\t<0> = 83.52 and (\ff<(CG)Symbol>m\ff<(CG)Times>/\ff<(CG)Symbol>s\ff<(CG)Times>)\SP2\Sp,\ll<1> calculated as
\t<-1>Reason for termination: gradient size is less than gradient tolerance.\t<0>
(83.5173 1.67099)\b\ff<(CG)Times>
Some statust information is printed as the optimization proceeds. You can turn this off by supplying the keyword argument \B\ff<(CG)Courier>:verbose\b\ff<(CG)Times> with value \B\ff<(CG)Courier>NIL\b\ff<(CG)Times>.
\piYou might want to check that the gradient of the function is indeed close to zero. If you do not have a closed form expression for the gradient you can use \B\ff<(CG)Courier>numgrad\b\ff<(CG)Times> to calculate a numerical approximation. For our example,\pn
The elements of the gradient are indeed close to zero. You can also compute the second derivative, or Hessian, matrix using \B\ff<(CG)Courier>numhess\b\ff<(CG)Times>. Approximate standard errors of the maximum likelihood estimates are given by the square roots of the diagonal entries of the inverse of the negative Hessian matrix at the maximum:
\piInstead of calculating the maximum using \B\ff<(CG)Courier>newtonmax\b\ff<(CG)Times> and then calculating the derivatives separately you can have \B\ff<(CG)Courier>newtonmax\b\ff<(CG)Times> return a list of the location of the maximum, the optimal function value, the gradient and the Hessian by supplying the keyword argument \B\ff<(CG)Courier>:return-derivs\b\ff<(CG)Times> as \B\ff<(CG)Courier>T\b\ff<(CG)Times>.\SP28\Sp
Newton's method assumes a function is twice continuously differentiable. If your function is not smooth or you are having trouble with \B\ff<(CG)Courier>newtonmax\b\ff<(CG)Times> for some other reason you might try a second maximization function, \B\ff<(CG)Courier>nelmeadmax\b\ff<(CG)Times>. \B\ff<(CG)Courier>nelmeadmax\b\ff<(CG)Times> takes a function and an initial guess and attempts to find the maximum using the Nelder-Mead simplex method as described in Press, Flannery, Teukolsky and Vetterling [14]. The initial guess can consist of a simplex, a list of \In\i+1 points for an \In\i-dimensional problem, or it can be a single point, represented by a list or vector of \In\i numbers. If you specify a single point you should also use the keyword argument \B\ff<(CG)Courier>:size\b\ff<(CG)Times> to specify as a list or vector of length \In\i the size in each dimension of the initial simplex. This should represent the size of an initial range in which the algorithm is to start its search for a maximum. We can use this method in our gamma example:\pn
It takes somewhat longer than Newton's method but it does reach the same result.
\BExercices\b
\po1. The data set used in this example consists of sets of measurements for ten aircraft. Data for five of the aircraft are contained in the variable failure-times in the file aircraft.lsp. The calculations of this section used the data for the second aircraft. Examine the data for the remaining four aircraft.
2. Use maximum likelihood estimation to fit a Weibull model to the data used in this section.
This section describes a set of tools for computing approximate posterior moments and marginal densities in XLISP-STAT. The definitions needed are in the file \B\ff<(CG)Courier>bayes.lsp\b\ff<(CG)Times> on the distribution disk. This file is not loaded automatically at start up; you should load it now, using the \B Load\b item on the \BFile\b menu or the \B\ff<(CG)Courier>load\b\ff<(CG)Times> command, to carry out the calculations in this section. The material in this section is somewhat more advanced as it assumes you are familiar with the basic concepts of Bayesian inference.
\piThe functions described in this section can be used to compute first and second order approximations to posterior moments and saddlepoint-like approximations to one dimensional marginal posterior densities. The approximations, based primarily on the results in [18, 19, 20], assume the posterior density is smooth and dominated by a single mode. The implementation is experimental and may change in a number of ways in the near future.
Let's start with a simple example, a data set used to study the relation between survival time (in weeks) of leukemia patients and white blood cell count recorded for the patients at their entry into the study [9, Example U]. The data consists of two groups of patients classified as AG positive and AG negative.\MThe data for the 17 AG positive patients, contained in the file \B\ff<(CG)Courier>leukemia.lsp\b\ff<(CG)Times> in the \B\ff<(CG)Courier>Data\b\ff<(CG)Times> folder on the distribution disk, can be entered as
\pnA high white blood cell count indicates a more serious stage of the disease and thus a lower chance of survival.
\piA model often used for this data assumes the survival times are exponentially distributed with a mean that is log linear in the logarithm of the white blood cell count. For convenience I will scale the white blood cell counts by 10000. That is, the mean survival time for a patient with white blood cell count log(\IWBC\i\SBi\Sb) is\pn
the log likelihood can be computed using the function
\B\ff<(CG)Courier>(defun llik-pos (theta)
\T (let* ((x transformed-wbc-pos)
\s (y times-pos)
\s (theta0 (select theta 0))
\s (theta1 (select theta 1))
\s (t1x (* theta1 x)))
\T (- (sum t1x)
\s(*\M(length x)\M(log theta0))
\s(/ (sum (* y (exp t1x)))
\s (theta0))))\b\ff<(CG)Times>
\piI will look at this problem using a vague, improper prior distribution that is constant over the range \ff<(CG)Symbol>q\SB\ff<(CG)Times>i\Sb > 0; thus the log posterior density is equal to the log likelihood constructed above, up to an additive constant. The first step is to construct a Bayes model object using the function \B\ff<(CG)Courier>bayes-model\b\ff<(CG)Times>. This function takes a function for computing the log posterior density and an initial guess for the posterior mode, computes the posterior mode by an iterative method starting with the initial guess, prints a short summary of the information in the posterior distribution, and returns a model object. We can use the function \B\ff<(CG)Courier>llik-pos\b\ff<(CG)Times> to compute the log posterior density, so all we need is an initial estimate for the posterior mode. Since the model we are using assumes a linear relationship between the logarithm of the mean survival time and the transformed \IWBC\i variable, a linear regression of the logarithm of the survival times on \B\ff<(CG)Courier>transformed-wbc-pos\b\ff<(CG)Times> should provide reasonable estimates. The linear regression gives
\pnso reasonable initial estimates of the mode are \ff<(CG)Symbol>q\SB\ff<(CG)Times>0\Sb = exp(3.5) and \ff<(CG)Symbol>q\SB\ff<(CG)Times>\ll<1>1\Sb = 0.8. Now we can use these estimates in the \B\ff<(CG)Courier>bayes-model\b\ff<(CG)Times> function:
It is possible to suppress the summary information by supplying \B\ff<(CG)Courier>NIL\b\ff<(CG)Times> as the value of the \B\ff<(CG)Courier>:print\b\ff<(CG)Times> keyword argument.
\piThe summary printed by \B\ff<(CG)Courier>bayes-model\b\ff<(CG)Times> gives first order approximations to the posterior means and standard deviations of the parameters. That is, the means are approximated by the elements of the posterior mode and the standard deviations by the square roots of the diagonal elements of the inverse of the negative Hessian matrix of the log posterior at the mode. These approximations can also be obtained from the model by sending it the \B\ff<(CG)Courier>:1stmoments\b\ff<(CG)Times> message:\pn
The result is a list of two lists, the means and the standard deviations. A slower but more accurate second order approximation is available as well. It can be obtained using the message \B\ff<(CG)Courier>:moments\b\ff<(CG)Times>:
Both these messages allow you to compute moments for individual parameters or groups of parameters by specifying an individual parameter index or a list of indices:
\piThe first and second order approximations to the moments of \ff<(CG)Symbol>q\SB\ff<(CG)Times>0\Sb are somewhat different; in particular the mean appears to be somewhat larger than the mode. This suggests that the posterior distribution of this parameter is skewed to the right. We can confirm this by looking at a plot of the approximate marginal posterior density. The message \B\ff<(CG)Courier>:margin1\b\ff<(CG)Times> takes a parameter index, and a sequence of points at which to evaluate the density and returns as its value a list of the supplied sequence and the approximate density values at these points. This list can be given to plot-lines to produce a plot of the marginal density:
The result is shown in Figure 20 and does indeed show some skewness to the right.
In additioin to examining individual parameters it is also possible to look at the posterior distribution of smooth functions of the parameters.\SP29\Sp For example, you might want to ask what information the data contains about the probability of a patient with \IWBC\i = 50000 surviving a year or more. The probability is given by\ll<-1>
\jcFigure 20:\MPlot of marginal posterior density for \ff<(CG)Symbol>q\SB\ff<(CG)Times>0\Sb.\jf
\ll<1>This function can be given to the \B\ff<(CG)Courier>:1stmoments\b\ff<(CG)Times>, \B\ff<(CG)Courier>:moments\b\ff<(CG)Times> and \B\ff<(CG)Courier>:marginal\b\ff<(CG)Times> methods to approximate the posterior moments and marginal posterior density of this function. For the moments the results are\ll<-2>
is shown in Figure 21. Based on this picture the data suggests that this survival probability is almost certainly below 0.5, but it is difficult to make a more precise statement than that.
\piThe functions described in this section are based on the optimization code described in the previous section. By default derivatives are computed numerically. If you can compute derivatives yourself you can have your log posterior function return a list of the function value and the gradient or a list of the function value, the gradient and the Hessian matrix.\pn
\BExercises\b
\po1. To be able to think about prior distributions for the two parameters in this problem we need to try to understand what the parameters represent. The parameter \ff<(CG)Symbol>q\SB\ff<(CG)Times>0\Sb is fairly easy to understand: it is the mean survival time for patients with \IWBC\i = 10000. The parameter \ff<(CG)Symbol>q\SB\ff<(CG)Times>1\Sb is a little more difficult to think abour. It represents the approximate percent difference in mean survival time for patients with \IWBC\i differing by one percent. Because of the minus sign in the mean relationship, and the expected inverse relation between \IWBC\i and mean survival time, \ff<(CG)Symbol>q\SB\ff<(CG)Times>1 \Sb is expected to be positive.
\sConsider an informative prior distribution that assumes the two parameters \Ia priori\i independent, takes log(\ff<(CG)Symbol>q\SB\ff<(CG)Times>0\Sb) to be normally distributed with mean log(52) and standard deviation log(2), and \ff<(CG)Symbol>q\SB\ff<(CG)Times>1\Sb to be exponentially distributed with mean \ff<(CG)Symbol>m\ff<(CG)Times> = 5. This prior is designed to represent an opinion that mean survival time at \IWBC\i = 10000 should be around one year, but that guess could easily be off by a factor of two either way. The percentage change in the mean for a one percent change in \IWBC\i should be on the order of one to ten or so. Examine the posterior distribution for this prior and compare your results to the results for the vague prior used above.
2. Construct and examine a posterior distribution for the parameters of the gamma model based on the aircraft data of Section 12.\pn
\B\fs<14>References\b\fs<12>
\T\po[1] Bates, D.M. and\MWatts, D.G., (1988), \INonlinear Regression Ananlysis and its Applications\i, New York: Wiley.
\T[2] Becker, Richard A., and Chambers, John M., (1984), \IS: An Interactive Environment for data Analysis and Graphics\i, Belmont, Ca: Wadsworth.
\T[3] Becker, Richard A., Chambers, John M., and Wilks, Allan R., (1988), \IThe New S Language: A Programming Environment for Data Analysis and Graphics\i,'' Pacific Grove, Ca: Wadsworth.
\T[4] Becker, Richard A., and William S. Cleveland, (1987), ``Brushing scatterplots'', \ITechnometrics\i, vol. 29, pp. 127-142.
\T[6] Betz, David, (1988), ``XLISP: An experimental object-oriented programming language,'' Reference manual for XLISP Version 2.0.
\T[7] Chaloner, Kathryn, and Brant, Rollin, (1988) ``A Bayesian approach to outlier detection and residual analysis'', \IBiometrika\i, vol. 75, pp. 651-660.
\T[8] Cleveland, W. S. and McGill, M. E., (1988) \IDynamic Graphics for Statistics\i, Belmont, Ca.: Wadsworth.
\T[9] Cox, D.R. and Snell, E.J. (1981) \IApplied Statistics:\MPrinciples and Examples\i, London: Chapman and Hall.
[10] Dennis, J.E> and Schnabel, R.B. (1983)\M\INumerical Methods for Unconstrained Optimization and Nonlinear Equations\i, Englewood Cliffs, N.J>:\MPrentice-Hall.
[11] Devore, J. and Peck, R., (1986), \IStatistics, the Exploration and Analysis of Data\i, St. Paul, Mn: West Publishing Co.
[12] McDonald, J. A., (1982), "Interactive Graphics for Data Analysis," unpublished Ph. D. thesis, Stanford University, Department of Statistics.
[13] Oehlert, Gary W., (1987), ``MacAnova User's Guide,'' Technical Report 493, School of Statistics, University of Minnesota.
[14] Press, Flannery, Teukolsky and Vetterling, (1988), \INumerical Recipes in C\i, Cambridge: Cambridge University Press.
[15] Steele, Guy L., (1984), \ICommon Lisp: The Language\i, Digital Press.
[16] Stuetzle, W., (1987), "Plot windows," \IJ. Amer. Statist. Assoc\i., vol. 82, pp. 466 - 475.
[17] Tierney, Luke, (1990) \ILISP-STAT An Object-Oriented Environment for Statistical Computing and Dynamic Graphics\i. New York: Wiley.
[18] Tierney, L. and Kadane, J.B., (1986), "Accurate approximations for posterior moments and marginal densities."\M\IJ. Amer. Statist. Assoc.,\i vol. 81, pp.82-86.
[19] Tierney, Luke, Kass, Robert, E. and Kadane, Joseph B., (1989), "Fully exponential Laplace approximations to expectations and variances of nonpositive functions." \IJ. Amer. Statist. Assoc.\i to appear.
[20] Tierney, Luke, Kass, Robert, E. and Kadane, Joseph B., (1989), "Approximate marginal densities for nonlinear functions,"\M\IBiometrika\i, to appear.
[21] Weisberg, Sanford, (1982), "MULTREG Users Manual," Technical Report 298, School of Statistics, University of Minnesota.
[22] Winston, Patrick H. and Berthold K. P. Horn, (1988), \ILISP\i, 3rd Ed., New York: Addison-Wesley.\pn
\SP\pi23\SpThe discussion in this section only scratches the surface of what you can do with functions in the XLISP language. To see more examples you can look at the files that are loaded when XLISP-STAT starts up. For more information on options of function definition, macros, etc. see the XLISP documentation and the books on Lisp mentioned in the\pn references.
\SP\pi24\Spmapcar\b\ff<(CG)Times> can be given several lists after the function. The function must take as many arguments as there are lists. \B\ff<(CG)Courier>mapcar\b\ff<(CG)Times> will apply the function using the first element of each list, then using the second element, and so on, until one of the lists is exhausted, and return a list of the results.\pn
\SP\pi26\SpRecall that the expression \B\ff<(CG)Courier>#'f1\b\ff<(CG)Times> is short for \B\ff<(CG)Courier>(function f1)\b\ff<(CG)Times> and is used for obtaining the function definition\pn associated with the symbol \B\ff<(CG)Courier>f1\b\ff<(CG)Times>.
\SP\pi27\b\SpThe maximizing value for \ff<(CG)Symbol>m\ff<(CG)Times> is always the sample mean. We could take advantage of this fact and reduce the problem to a one dimensional maximization, but it is simpler to just maximize over both parameters.\ff<(CG)Courier>\pn
\SP\pi28\SpThe function \B\ff<(CG)Courier>newtonmax\b\ff<(CG)Times> ordinarily uses numerical derivatives in its computations. Occasionally this may not be accurate enough or may take too long. If you have an expression for computing the gradient or the Hessian then you can use these by having your function return a list of the function value and the gradient, or a list of the function value, the gradient and the Hessian matrix, instead of just returning the function value.\pn
\n\ff<(CG)Times>\fs<12>\jc\ll<-2>\ps<100>\t<0>\FpFigure 21: Plot of marginal posterior density of the one year survival probability of a patient with \IWBC\i = 50000.\.