Introduction to Data Analysis with Mathematica

1. Importing and Exporting Data

Preliminaries

Strings

A string is a sequence of characters inside quotation marks. Because Mathematica does not normally interpret the contents of a string as Mathematica commands, a string is a useful way to represent many kinds of data. For example, file and directory names are represented as strings.

There are several special characters that cannot be directly entered in a string, except by using escape sequences. String escape sequences consist of a backslash followed by a single letter.

Mathematica treats a newline between quotation marks as a space.

"Sentence one.
	Sentence two."

To put a proper newline character in a string, we use the escape sequence \n.

"This is line one.\nAnd this is line two."

The input form of the string displays the escape sequence used to enter the newline character.

InputForm[%]

To put a quotation mark into a string (instead of using it to denote the beginning or end of a string), we use the escape sequence \".

"We enter quotation marks using an \"escape sequence\"."

Again the input form reveals any escape sequences contained in a string.

InputForm[%]

Because the backslash character is used to mark the beginning of an escape sequence, to enter a literal backslash character into a string we must use the escape sequence \\. This is important when using a file system such as Microsoft Windows that uses a backslash as a directory separator. A typical directory name on such a system looks like the following. The output form of the string displays a single backslash for each double backslash contained in the input.

"C:\\data"

Here is the input form of the string.

InputForm[%]

Many strings can be converted into valid Mathematica input. Mathematica does not evaluate the contents of the following string.

"1 + 1"

The function ToExpression converts a string into Mathematica input, then evaluates it, if possible.

ToExpression[%]

We will see in a later chapter that it is occasionally useful to import numerical data into Mathematica as a list of strings, then convert the strings into numbers using ToExpression.

Numbers

Mathematica has different rules for working with exact and approximate quantities.

An approximate quantity can be identified by the presence of a decimal point; therefore 17.0 and 0.71 are approximate numbers. The integer 17 and the ratio 71/100, on the other hand, do not contain a decimal point and are considered exact numbers. In the following discussion, we refer to an approximate (non-complex) number containing a decimal point as a real number, an integer without a decimal point as an integer, and a ratio of two exact integers as a rational number. To Mathematica 7.0 is an approximate real number, and 7 is an exact integer.

When doing arithmetic with exact numbers, Mathematica returns an exact number. For example, the following input contains only exact integers and rational numbers, so the result is an exact number.

2 + (((3×(5/12)))^7)

This occurs even when there is a precise decimal equivalent.

(2/3) + (7/12)

If a calculation involves even one approximate number (a number that contains a decimal point), however, the result will be an approximate number because the uncertainty associated with the approximate number propagates through the whole calculation.

(2/3) + (7.0/12)

I denotes the imaginary unit (Sqrt[-1]).

Sqrt[-256]

Mathematica knows the standard functions for describing and manipulating complex numbers. Re[z] returns the real part of z, Im[z] returns the imaginary part, and so forth.

Here is a complex number called num.

num = 1.23 + 4.56 I

Here is a list of the real and imaginary parts of num.

{Re[num],  Im[num]}

Here are the absolute value and approximate argument (or phase), in radians, of num.

{Abs[num],  Arg[num]}

To convert from an exact number to an approximate number, we use the function N. N[expr] returns a numerical approximation to expr.

N[(Sqrt[27]) - (Sqrt[13])]

N[expr, n] does computations to at most n significant digits. Here is a 50-digit approximation to the difference between (22/7) and π.

N[22/7 - π, 50]

The presence of the decimal point in the following result indicates the change from exact integer to approximate real number.

N[2]

There are many functions that convert an approximate number into an exact integer. Round[x] returns the integer closest to x, Floor[x] returns the greatest integer less than x, and Ceiling[x] returns the least integer greater than x.

{Round[3.3], Floor[3.3], Ceiling[3.3]}

The function Rationalize converts numbers into exact rational numbers. Rationalize[x,dx] returns a rational number equal to x within a tolerance of dx.

Here are some successively better rational approximations to π.

Table[Rationalize[π, (10^(-(2^n)))], {n, 1, 4}]

We can find the type of a number using the function Head.

{Head[2],  Head[3.2], Head[3 + 4  I], Head[2/3]}

Mathematica can work with numbers with any number of digits, and in the case of approximate numbers it maintains information about how many of the digits are significant. The number of significant digits in a number x is called the precision of x, and Mathematica knows the rules of numerical analysis for dealing with numbers with different amounts of precision.

Precision[x] returns the number of significant digits in the number x.

Precision[N[π, 40]]

Accuracy returns the number of digits to the right of the decimal point. Precision is a measure of the relative error of a number, and accuracy is a measure of the absolute error.

The following number has 24 significant digits.

Precision[1234567890.123456789012345]

It has 15 digits to the right of the decimal point.

Accuracy[1234567890.123456789012345]

Machine-precision numbers are numbers that can be calculated using a computer’s hardware: most computers can directly handle numbers with up to 16 significant digits. (The number of digits that the hardware can handle is stored in the parameter $MachinePrecision.) For efficiency reasons, Mathematica does any calculation that contains even one machine-precision number to machine precision.

The number 2.3 in the following calculation is a machine-precision number (while the other number is not), so the answer is a machine-precision number.

2.3 + 1.234567890123456789012345

Precision[%]

By default, N returns a machine-precision number.

N[(Sqrt[19])]

Precision[%]

If a number has more than 16 (or the value of $MachinePrecision, if different) significant digits, it is called an arbitrary-precision number.

Here we ask for an approximation to (Sqrt[19]) with 50 significant digits.

N[(Sqrt[19]),  50]

The precision of the result is 50, making the approximation an arbitrary-precision number.

Precision[%]

Using SetPrecision we can artificially set the precision of a number to an arbitrary number of digits.

SetPrecision[ N[(Sqrt[19]), 4],  5]

Precision[%]

When using arbitrary-precision numbers, Mathematica does not generate more precision in an answer than is justified by the calculation. For example, adding two values with different numbers of significant digits returns a number as precise as the less precise value.

\(x\_1 = 54.23232323232312312624874590149643; 

    x\_2 = 34323.98129712872939137913; \)

Precision[(Subscript[x, 1])]

Precision[(Subscript[x, 2])]

Precision[(Subscript[x, 1]) + (Subscript[x, 2])]

The precision of an exact number is Infinity.

We can approximate an exact number to as many digits as desired.

N[Sin[1], 250]

However, we cannot approximate the following machine-precision number to more than machine precision.

N[Sin[1.0], 250]

An additional way to enter a number with a known precision is in the form nnnn`p, where nnnn is the number and p is the number of digits of precision. For example, here is an approximation to π with four digits of precision.

3.1416`4

By default Mathematica prints only the correct digits, without the `p. To explicitly see all of the information Mathematica has for a number, we look at the number’s input form.

InputForm[%]

As with arbitrary-precision arithmetic, Mathematica keeps track of the number of significant digits in a result. Here we use several imprecise numbers in the calculation of the volume of a cylinder, and Mathematica returns an answer with the correct number of significant digits.

3.1416`4*(((11.1111`5))^2) (15.253545`6)

Where Mathematica Looks for Files

Mathematica maintains a current working directory, which is the first directory in which Mathematica looks for files. We can find the current working directory by entering the command Directory[].

Directory[]

The current directory will appear somewhat different on different computers because the delimiters used to separate directories differ from system to system. This notebook was created on a Macintosh, which uses “:” to separate directories in a pathname. The delimiter is “/”  in Unix and “\” in Microsoft Windows; the pathname separator for a particular computer can be found by looking at the value of the global parameter $PathnameSeparator.

If we refer to a file is that not in the current working directory, Mathematica looks for the file in the directories listed in the variable $Path. This list typically includes the Mathematica application directory and the StandardPackages subdirectory within it. Note that the value of $Path may be different on different computers.

$Path

Mathematica searches the directories listed in $Path in the order they are listed.

The value of $Path is a list. For this training course, we assume the data files are in a directory called data, located at the top level of the file system. The following line of code adds this directory to the path on a Macintosh system.

AppendTo[$Path, "HardDisk:data"]

On a Unix system, we enter AppendTo[$Path,"/data"]; and on a Microsoft Windows system, we enter AppendTo[$Path,"c:\\data"]. (The double backslash is explained earlier in the section “Strings”.)

Because $Path is a list, we treat it just like any other list in Mathematica, using list commands (such as AppendTo) to change its contents.

Here we use SetDirectory to set the current directory to be the directory in which our data files are stored.

SetDirectory[Last[$Path]]

Referring to Directories and Files

Mathematica contains many commands that enable us to interact with the local file system. We can view the contents of directories and files, create new directories and files, read the contents of files into Mathematica data structures, and much more.

FileNames[ ] lists all the files in the current working directory. Here is a list of every file in the current directory (HardDisk:data).

FileNames[]

FileNames["form"] lists all files in the current working directory whose names match the string pattern form. Here we use the string pattern “*m*” to find all file names in the current directory that contain the letter m.

FileNames["*m*"]

Here we set the current directory to the NumericalMath add-on package directory. (When using a notebook interface to Mathematica, we can pull down the Input menu and choose Get File Path, which presents a standard file dialog box. Using the dialog box, we can navigate directories until we find a file in the desired directory. Mathematica pastes the complete pathname to the file at the current text insertion point. To change the file pathname to a directory name, we simply delete the file name from the complete pathname.)

SetDirectory[
      "HardDisk:Mathematica 3.0 Files:AddOns:StandardPackages:NumericalMath\
"]

Here is a list of all the filenames in the NumericalMath directory.

FileNames[]

FindList["file", "string"] gives a list of all records in file that contain string. Here we find every line in the file TrigFit.m that contains the word “Author”.

FindList["TrigFit.m", "Author"]

Mathematica keeps a list of current directories. ResetDirectory steps backward through this list.

ResetDirectory[]

Exercises: Preliminaries

Make sure that the data file directory is on the file search path $Path. (The directory is called data, located on the root level of the file system.)

If the data directory is not listed in $Path, we append it to the list by using AppendTo.

AppendTo[$Path, "HardDisk:data"]

The data directory is now at the end of the search path.

$Path

Set the current directory to the data directory.

Here we set the directory using SetDirectory and the full pathname to the data directory.

SetDirectory["HardDisk:data"]

Alternatively, if the data directory is the last directory listed in $Path, we can enter the following.

SetDirectory[Last[$Path]]

List all the files in the data directory.

FileNames[] returns a list of all the files in the current directory.

FileNames[]

Importing Data

Viewing Files

The command !!file prints the contents of a file into the current notebook. Here are the contents of the file called sample.dat inside the data directory.

!!sample.dat

To access a file without adding its path to $Path, we use the full pathname. Here is the full pathname for the file sample.dat.

!!HardDisk:data:sample.dat

Recall that pathname directory separators are different for different filesystems. We enter the equivalent pathname on a Microsoft Windows system as “C:\\data\\sample.dat”, and on a Unix system as “/data/sample.dat”.

Reading Lists of Data

The basic command used to import data from a file is ReadList. By giving suitable values to the options that ReadList accepts, we can read data stored in almost any format.

ReadList["file"] returns a list of all the expressions in file. Here we use ReadList to read the contents of the file sample.dat.

ReadList["sample.dat"]

We tell Mathematica what type of data to import by specifying a second argument to ReadList. ReadList["file", type] reads objects of the specified type from file, until the end of file is reached. ReadList["file",{(Subscript[type, 1]),(Subscript[type, 2]),…}] reads objects of the given sequence of types, until the end of file is reached.

Number, Word, and String are the most common types of data to read; other types include Byte, Character, Real, Record, and Expression. The important thing to remember is that once the data are read into a Mathematica list, we manipulate the list using the same functions as any other Mathematica data structure.

Here we read the contents of the file sample.dat as a list of numbers.

data = ReadList["sample.dat", Number]

The first element of the list is an integer.

Head[First[data]]

The last element is a real number.

Head[Last[data]]

Here we read the same file as a list of words.

data = ReadList["sample.dat", Word]

Each element of the new list is a string.

Head[First[data]]

It is good practice to clear the value of a variable when it is no longer needed.

Clear[data]

A third argument to ReadList specifies how many objects to read. Here we read the first three elements of the file sample.dat.

ReadList["sample.dat", Number, 3]

Mathematica also reads data files in which expressions are separated by tab characters.

!!tabs.dat

ReadList["tabs.dat", Number]

ReadList interprets numbers of the form 123E-2 as 123×(10^(-2)), using Fortran notation.

(Note that the form ReadList["file"], in which the type of data to be read is unspecified, treats each line in the file as a separate expression. Here we read the file tabs.dat without specifying that the elements are numbers.

ReadList["tabs.dat"]

Mathematica reads in each line as a separate expression, and returns the product of the numbers on each line.)

If a data file has a particular structure, it may be lost when using ReadList, since by default ReadList returns a simple list of the data points. The file lowtemps.dat contains two columns of data.

!!lowtemps.dat

When we read the data using ReadList, the column format is lost.

ReadList["lowtemps.dat", Number]

To preserve the structure, we alter the second argument to ReadList, in this case telling Mathematica to read in pairs of numbers.

ReadList["lowtemps.dat", {Number, Number}]

Here are the options that ReadList accepts.

Options[ReadList]

We use InputForm to reveal the input forms of the nonprinting characters in the option values. The valid record separator character is a newline character, whose input form is the escape sequence \n, and one of the valid word separators is a tab character, entered \t.

InputForm[%]

Alternatively, ReadList can be configured to create a sublist for each record. By default, ReadList considers a record to end with each newline character, and a word to end with each space. By setting the option RecordLists to True, we tell ReadList to create a new sublist each time a new record is encountered.

ReadList["lowtemps.dat", Number, RecordLists -> True]

Irregular arrays can also be read into Mathematica using ReadList, setting the option RecordLists to True. Again Mathematica collects data elements on the same line in a file into the same list.

!!primes.dat

ReadList["primes.dat", Number, RecordLists -> True] // 
      TableForm

Comma-Separated Data

By default, ReadList assumes that data are separated by spaces or tab characters, which is not always the case. In the file commas.dat, entries on the same line are separated by commas.

!!commas.dat

There are many techniques to read this data set into Mathematica. To read it as one flat list, we use the option WordSeparators to specify that words are separated by spaces or commas. Mathematica now treats commas in the data file as delimiters between words.

ReadList["commas.dat", Word, 
      WordSeparators -> {" ", ","}]

Now we can use Partition to recapture the lost structure. Here we recover the original structure of the data by partitioning the previous list into sublists of length four.

Partition[%, 4] // TableForm

Alternatively, we can collect elements on the same line into the same sublist by setting the option RecordLists to True.

commadata = 
      ReadList["commas.dat", Word, 
          WordSeparators -> {" ", ","}, RecordLists -> True]
           // TableForm

Because we read in the data as words, the numbers contained in the data file are represented as strings.

InputForm[First[commadata]]

Here we use ToExpression to convert the strings containing numeric data into numbers.

commadata = ToExpression[commadata]

The numbers in the data file are now represented as numbers.

InputForm[First[commadata]]

Here we clear the value of commadata.

Clear[commadata]

Files Containing Numerical and Textual Data

Many data files contain both numerical and textual data. The following file contains a list of geological periods and their respective durations.

!!periods.dat

The file contains a list of records, each containing a word and a number. We must supply Mathematica with this information in order to read in the data correctly. In this case, we simply tell Mathematica to read the data into the appropriate list structure.

ReadList["periods.dat", {Word, Number}]

Each pair of elements is now contained in a list containing a word and a number. We see that Mathematica internally stores data read in as words in the form of strings.

InputForm[%]

Here is a similar data file.

!!rankbypop.dat

Just as before, we want to read in the file as a list of records of the form {state, population-rank}. However, if we use the same command as for the previous data file, Mathematica prints an error message.

ReadList["rankbypop.dat", {Word, Number}]

The problem is caused by the space between the words “New” and “York”. Mathematica reads the string “New York” as two data items, and prints a message that “York” cannot be interpreted as a number.

In this case, we can avoid this problem by specifying that words are separated by two spaces, using the option WordSeparators.

\(ReadList["rankbypop.dat", {Word, Number}, 
	
      WordSeparators -> "  "\  (*two\ spaces*) ]\)

In other cases, more sophisticated file-reading techniques may be required.

Streams

Intstead of reading an entire file at once, we can open a channel to a file, called a stream. When we open a stream to a file, Mathematica holds open a connection to the file, and allows us more flexibility in reading from and writing to the file. OpenRead["file"] opens a stream to read data from the file named file, and returns an InputStream object. We name the stream in order to make it easier to refer to the stream.

practicestream = OpenRead["sample.dat"]

(The number inside the InputStream object is a unique index that Mathematica generates to distinguish streams.)

Rather than reading the entire file into a list, here we read only one number from the file.

Read[practicestream, Number]

Here we create a table of values read separately from the file. Mathematica returns the symbol EndOfFile if we try to read past the end of a data file.

Table[Read[practicestream, Number], {4}]

Since Mathematica is holding the named file open, it is important to close the file when we are finished reading from it using Close.

Close[practicestream]

Streams allow a great degree of control when dealing with files. For example, we can move within the file without reading the entire file into Mathematica. Read[stream, type] reads one expression of the specified type from a stream, and Read[stream, {type1, type2, …}] reads a sequence of objects of the specified types. Skip[stream, type, n] skips n objects of the specified type.

When used with a stream, ReadList reads in the remaining contents of a file, from the current position in the file to the end of the file.

Here is a file with some lines of header information before the data.

!!planets.dat

We cannot use ReadList alone to read this information in because the structure of the data is not uniform. What we need to do is skip the first three lines, and then apply ReadList to the remaining lines of the file. Here we open a stream to the file containing the planetary data.

planetstr = OpenRead["planets.dat"]

Then we use the Skip command to skip past three objects of type String. In a file or stream, an object of type String is an entire line of text.

Skip[planetstr, String, 3]

StreamPosition returns the current position in a stream, typically the number of bytes from the beginning of the stream. (We can move the current position using SetStreamPosition. To move the current position to the beginning of a stream, for example, we type SetStreamPosition[stream,0].)

StreamPosition[planetstr]

Once we have skipped the first three lines, the data follow a structured pattern which can be handled by ReadList. When we apply ReadList to an open stream, Mathematica reads only the contents of the file after the current position in the stream.

ReadList[planetstr, {Word, Number}]

Note that the file remains open until we close the stream, preventing other applications from accessing it.

Close[planetstr]

Exercises: Importing Data

Read in the file primes.dat first as a simple list of numbers, then as a list of pairs of numbers.

Here is the data as a simple list of numbers.

ReadList["primes.dat", Number]

Here we read in the data as a list of pairs of numbers.

ReadList["primes.dat", {Number, Number}]

Notice that there is an odd number of data points in the file, so the last list of pairs contains the symbol EndOfFile, which indicates that the end of the file was reached, so a second number in the pair could not be read.

Exporting Data

Writing to Files

When writing data to files using the Mathematica programming language, we use similar structures to those used for reading data from files.

Put and PutAppend write expressions or sequences of expressions to a file.

Put[Table[Random[], {5}], "rand5.dat"]

Here are the contents of the file rand5.dat.

!! rand5.dat

The following is an abbreviation for Put.

Table[Random[], {5}] >> rand5.dat

We can write data incrementally to an output stream, instead of writing an entire data structure to a file at once. Write[stream, (Subscript[expr, 1]), (Subscript[expr, 2]), …] writes the expressions (Subscript[expr, i]) in sequence, followed by a newline character, to the specified stream; WriteString[stream, (Subscript[expr, 1]), (Subscript[expr, 2]), …] converts the (Subscript[expr, i]) into strings, but does not append the newline character.

OpenWrite["file"] creates and opens a file for writing, returning an OutputStream object. This will destroy any data in an existing file. To append data to the end of an existing file, we use OpenAppend["file"].

Here we use OpenWrite to create a file called newfile. (OpenWrite deletes the contents of an existing file.)

newstream = OpenWrite["newfile"]

Here we write a list to the stream.

Write[newstream, {{2, 1, 2}, {1, 2, 1}}]

Here we close the stream.

Close[newstream];

Here are the current contents of newfile.

!! newfile

To append to the contents of a file, we use OpenAppend.

addtofile = OpenAppend["newfile"]

Here we write another expression to newfile.

Write[addtofile, a x^2 + b x + c == 0]

Here we write some text to the file.

WriteString[addtofile, 
	"\n\nThe preceding lines were written\n", 
	"using Write, but these were added\n",
	"with WriteString."]

When finished writing to newfile, we close the stream.

Close[addtofile]

These are now the contents of newfile.

!!newfile

Writing Arrays to Files

Passing an array structure though Write or WriteString produces formatted output, including the curly brackets that Mathematica uses to denote lists, as the following example demonstrates. Here is a data set.

somedata
       = {{3.50, 77, 83, 90, 98}, {3.75, 82, 89, 97}, {4.00, 88, 95}, {4.25, 
          93, 100}, {4.50, 99}};

Here we write the data set directly to a file.

arraytest = OpenWrite["array.txt"]

Write[arraytest, somedata]

Close[arraytest];

The contents of the file follow Mathematica syntax conventions.

!! array.txt

To write a formatted array, one possibility is to format the array using TableForm, then write the formatted array to a file.

WriteString["array1.txt", 
      OutputForm[TableForm[somedata, TableSpacing -> {0, 3}]]]

Here are the contents of our file.

!! array1.txt

Here we read the data back into Mathematica.

ReadList["array1.txt", Number, RecordLists -> True]

Similarly, we can use the following function definition. (We begin the name writeArray with a lower case letter to distinguish it from built-in functions.) The nested Do loops step through the array row by row, then column by column. Each element of a row is written followed by a space. Using WriteString ensures that there are no newline characters inserted. Between rows Write is used to insert a newline character.

\(\(writeArray[streamname_, data_List] := 
	
      Do[
		Do[
			
          WriteString[streamname, data[\([i, j]\)], " "], {j, 1, 
            Length[data[\([i]\)]]}]; 
        Write[streamname], {i, 1, Length[data]}]; \)\)

(An equivalent function can be written using pure functions and the Scan command.)

Here is a test of the function.

arraytest2 = OpenWrite["array2.txt"]

writeArray[arraytest2, somedata]

Close[arraytest2];

Here are the contents of the file written by writeArray.

!! array2.txt

Exercises: Exporting Data

Create a new file and write a few lines of text to it.

Here we create a file called export, and open it for writing. Mathematica returns an open stream to the file.

myfile = OpenWrite["export"]

Here we write some text to the file.

WriteString[myfile, "Here are the first ten prime numbers:\n"]

Note that WriteString converts any expression into string form before writing it to a file. Here we write the first ten prime numbers, with a space after each number, to a file using a Do loop.

Do[WriteString[myfile, Prime[i], " "], {i, 1, 10}]

Here we add two newline characters to the end of the file.

WriteString[myfile, "\n\n"]

When we are finished writing to the file, we close it.

Close[myfile]

Here are the contents of the file.

!!export

Create a 10 × 10 array of random numbers of any type or size. Using the writeArray function defined in the text, append the array to the file export.

Here is a 10 × 10 array of random integers between one and ten.

ran = Table[Random[Integer, {1, 10}], {10}, {10}];

Here is the definition of writeArray from the discussion.

\(Clear[writeArray]; 

    writeArray[streamname_, data_List] := 
      Do[Do[WriteString[streamname, data[\([i, j]\)], " "], {j, 1, 
            Length[data[\([i]\)]]}]; 
        Write[streamname], {i, 1, Length[data]}]; \)

OpenAppend open an output stream to a file.

myfile = OpenAppend["export"];

Here we write the random array to the stream.

writeArray[myfile, ran]

Here we close the stream.

Close[myfile];

These are now the contents of the file.

!!export

Clear[ran]

Executing Files

If a file contains Mathematica commands, we often wish to execute the contents of the file. When we have valid Mathematica expressions in a file, Get["file"] reads in the file, evaluates each expression in it, and returns the result of the last expression. <<file is an abbreviation for Get["file"]. Needs["context`"] is used to read in Mathematica packages, and loads a package only if the context defined in the package does not already exist.

An example is the file execute, which contains a matrix definition. Here are the contents of execute.

!!execute

Get evaluates the contents of execute, and returns the result of the evaluation.

Get["execute"]

Here we see the value of newdata has been assigned.

newdata

Here we calculate the maximum value contained in newdata.

Max[newdata]

Here we execute the contents of a larger file called eye.dat, which defines a list called eye to be a list of gray values taken from a digitized photograph.

Get["eye.dat"];

Here is a density plot of the data.

ListDensityPlot[eye, Mesh -> False, 
      AspectRatio -> Automatic];

To save the definition of a symbol into a Mathematica-executable file, we use Save. Here is a list of data points.

circledata = Table[{Cos[t], Sin[t]}, {t, 0, 1.99  π, 0.5}]

Here we save the definition of circledata into a file.

Save["circle.dat", circledata]

Here are the contents of the file circle.dat.

! ! circle.dat

The contents of the file circle.dat are now executable using Get. First we clear the value of the variable circledata.

Clear[circledata]

We now execute the contents of circle.dat, which restores the definition of circledata.

Get["circle.dat"];

The results are now usable in further calculations.

ListPlot[circledata, PlotStyle -> PointSize[0.03], AspectRatio -> 1];

The function DumpSave works indentically to Save, except that the files are saved in a machine-specific binary format which loads faster than files created by Save. Here we save the data contained in the variable eye in Mathematica’s binary format. (By convention, files created by DumpSave have the filename extension .mx.)

DumpSave["eye.mx", eye];

To load the definition of eye contained in the binary file, we use Get. Note that files created by DumpSave are not compatible among different types of computers.

Clear[eye]

2. Lists

Creating and Interpreting Lists

The basic data structure in Mathematica is the list. The elements of a list can be any combination of numbers, variables, strings, functions, and other lists; enclosed within curly brackets {} and separated by commas.

Mathematica lists often have natural interpretations as vectors, matrices, or tensors, depending on how deeply they are nested. A one-dimensional list acts as a vector, a two-dimensional list as a matrix, and higher-dimensional lists as tensors.

One way to enter an arbitrary list is to type a sequence of elements inside curly brackets and separated by commas. The following one-dimensional list vec can be considered a vector with the three components a, b, and c.

vec = {a, b, c};

An easier way to create a vector of evenly spaced numbers is to use Range. Range[max] creates a list of integers from 1 through max; Range[min, max] creates a list of numbers from min through max in increments of 1; and Range[min, max, inc] returns a list of numbers ranging from min through max with increment inc.

vec2 = Range[2, 5, 1/2]

A more general function for creating lists is Table, which generates lists (nested to any depth) of elements that may depend on the given iterators. For example, to generate a list of values of (k^2) as k takes on values from 1 to 10, we evaluate the following.

Table[(k^2), {k, 1, 10}]

We create tables of more complicated functions in a similar way.

vec3 = 
      Table[N[Exp[π I t]], {t, 0, 2 π, π/4}]

Many mathematical and plotting functions are designed to operate on lists. Here is a plot of the real part of each element of vec3.

ListPlot[Re[vec3], PlotStyle -> PointSize[0.03]];

Fourier expects a list of values.

Fourier[vec3]

We generate two-dimensional tables by including two iterators.

rect = Table[(i/j), {i, 1, 5}, {j, 1, 2}]

It may be useful to use TableForm to display multidimensional lists.

TableForm[rect]

The options TableAlignments, TableSpacing, and TableHeadings allow us to customize the appearance of a table.

TableForm[rect, TableAlignments -> {Center, Center}, 
      TableSpacing -> {1, 10}, 
      TableHeadings 
        -> {
          None, {StyleForm["before", FontSize -> 9], 
            StyleForm["after", FontSize -> 9]}}]

How we interpret an array such as rect depends on the context in which it arises: in some cases it may be considered a 5 × 2 matrix, in others a list of five ordered pairs. Here we interpret rect as a matrix, and matrix-multiply it by its transpose using Dot (.).

rect.Transpose[rect] // MatrixForm

Here we interpret rect as a list of ordered pairs to be plotted.

ListPlot[rect, PlotStyle -> PointSize[0.02]];

The first argument to Table can be a conditional statement such as If or Which. An element of the following matrix is 1 if its indices i and j satisfy \(TraditionalForm
      \`|i - j| ≤ 2\), and is 0 otherwise.

MatrixForm[Table[If[Abs[i - j] ≤ 2, 1, 0], {i, 7}, {j, 7}]]

There is no need for Mathematica arrays to be rectangular. Here the upper limit of the inner loop j depends on the value of the outer loop i.

Table[i + j, {i, 1, 4}, {j, 1, i}]

The result is a triangular array.

TableForm[%]

There is no limit to the number of dimensions that can be contained in a list. Here we generate a rank-4 (four-dimensional) tensor by including four iterators in the Table command.

tensor = Table[i + j + k - l, {i, 2}, {j, 3}, {k, 4}, {l, 5}];

TableForm[tensor]

It is important to make a distinction between list operations and matrix operations, particularly with respect to matrix multiplication. The function Dot (which we enter as a period .) is used for matrix and vector multiplication. Here is the product of a matrix and a vector.

{{a, b}, {c, d}}.{x, y}

There is no difference in Mathematica between a row vector and a column vector, as matrix and vector functions automatically interpret vectors in the appropriate way. Here we multiply a matrix by a vector on both the left and the right; in both cases, we type the vector the same way.

{x, y}.{{a, b}, {c, d}}.{x, y} // Simplify

Here is the matrix product of two matrices.

{{1, 2}, {3, 4}}.{{a, b}, {c, d}}  // MatrixForm

Here we use ordinary multiplication *, and Mathematica performs termwise multiplication.

{{1, 2}, {3, 4}}*{{a, b}, {c, d}} // MatrixForm

As a general rule, if two lists have the same shape, Mathematica will perform termwise arithmetic, and return a list with the same shape. For example, the difference of two vectors of equal length is a vector of the differences of the elements.

{1, 2, 3} - {a, b, c}

Notice also that Mathematica will combine a list with a single (scalar) expression. If we add a constant to a vector, Mathematica returns a vector where the constant has been added to each element.

{a, b, c} + 2

Similarly, multiplying or adding a scalar to a matrix results in the expression being multiplied or added to each element of the matrix.

{{1, 2}, {3, 4}} 𝒞 //  MatrixForm

{{1, 2}, {3, 4}}  + 𝓀 //  MatrixForm

Mathematica produces a warning message if two lists cannot be combined because they have different lengths.

{1, 2, 3} + {a, b, c, d}

Exercises: Creating and Interpreting Lists

Create a list of the first ten prime numbers. (Prime[n] returns the (n^th) prime number.)

We use Table to generate a list of the (k^th) prime numbers as k goes from 1 to 10. When we omit a starting value for the iterator k, it is assumed to be 1.

Table[Prime[k], {k, 10}]

An alternative is to avoid the Table command by recognizing that the mathematical function Prime automatically maps itself over lists.

Prime[Range[10]]

Define a function called zeroVector that creates a vector containing n zeros, where n is an argument to zeroVector.

Here is one definition using Table. We can omit the name of the iterator if it does not occur in the expression given to Table.

zeroVector[n_Integer?Positive] := Table[0, {n}]

Here is a vector of ten zeros.

zeroVector[10]

We can define a zeroMatrix function in a similar way.

zeroMatrix[m_, n_] := Table[0, {m}, {n}]

Here is a 3 × 5 matrix of zeros.

zeroMatrix[3, 5] // MatrixForm

Generate a 5 × 5 matrix, in which each element is the row number plus the column number, or the number 5, whichever is greater.

We use Max to determine if 5 or i + j (the row number plus the column number) is greater.

Table[Max[{5, i + j}], {i, 5}, {j, 5}] // MatrixForm

Beginning with the matrix {{a, b}, {c, d}} and the vector {x, y}, create the vector {a x + b y, c x + d y}. Beginning with the same matrix and vector, create the matrix {{a x, b x}, {c y, d y}}.

This exercise illustrates the difference between list and matrix operations. Here are definitions of the matrix and vector.

exermat = {{a, b}, {c, d}}; exervec = {x, y};

For matrix-vector multiplication we use Dot.

exermat.exervec

To perform termwise multiplication we use Times.

exermat*exervec // MatrixForm

Compute the eigenvectors of a numeric approximation to the matrix

and see if they are normalized (that is, see if each eigenvector has unit length).

Here are the eigenvectors of the approximate matrix.

evecs=Eigenvectors[N[(({{1, 2, 3}, {3, 2, 2}, {7, 1, 2}}))]];TableForm[evecs]

The following function determines the length of a vector.

veclength[z_?VectorQ] := (Sqrt[z.z])

Here we use Table to test if each eigenvector has unit length.

Table[veclength[evecs[[i]]], {i, 3}]

Alternatively, we use the high-level programming function Map to map veclength over each vector in evecs.

Map[veclength, evecs]

Extracting Parts of Lists

Given a one-dimensional list, we can extract an element according to its position by using Part. In Mathematica, the first element of a list is at position 1, not 0.

Part[{a, b, c}, 1]

A compact, equivalent way to express the same command is to enclose the desired position between [[ and ]].

{a, b, c}[[1]]

In a multidimensional array, an element of a list can be another list. For example, here is a table of hurricane names, wind speeds (in knots), and pressures (in millibars), with each triple collected in a sublist.

hurr =
	{{"Allison",   65,  987},
	 {"Barry",     50,  995},
	 {"Chantal",   60,  992},
	 {"Dean",      40,  999},
	 {"Erin",      75,  975},
	 {"Felix",    115,  950},
	 {"Gabrielle", 60,  998},
	 {"Humberto",  90,  970},
	 {"Iris",      75,  970},
	 {"Jerry",     35, 1003}};

TableForm[hurr, TableSpacing -> {0, 3}]

To extract an entire sublist, or record, we specify the position of the record using Part (or the part-specification brackets [[ ]]).

Part[hurr, 2]

To extract a single element of a two-dimensional list, we refer to its coordinate, its row and column. For instance, to extract the wind speed of hurricane Erin, we first select Erin’s whole record, which is at position 5, and then select the second part of the record.

Part[hurr, 5, 2]

Here are two more ways to enter the same expression.

{hurr[[5, 2]], hurr[[5]][[2]]}

Similarly, we extract the pressure of the storm Gabrielle by extracting the sublist at position 7, then extracting the third part of Gabrielle’s record.

hurr[[7, 3]]

To extract elements from deeply nested lists, it may be necessary to specify several coordinates.

It is important to be aware that there is a difference between the part specification hurr[[1,2]] and hurr[[{1,2}]] (the latter has curly brackets inside the double square brackets). The former returns the element at position (1, 2) of the list hurr.

hurr[[1,2]]

The latter is equivalent to {hurr[[1]], hurr[[2]]}.

hurr[[{1, 2}]]

Note that [[{ }]] notation provides a quick way to manipulate the elements of a list. Here is an example of using this notation to rearrange the elements of a list.

plist = {a, b, c, d, e, f};

plist[[{1, 4, 2, 5, 3, 6}]]

Using Table and Part, we can define a function called sampleWithReplacement that takes a sample of a given size from a list of data.

sampleWithReplacement[L_List, n_Integer?Positive] := 
      Module[{ℓ = Length[L]}, 
        Table[L[[Random[Integer, {1, ℓ}]]], {n}]]

sampleWithReplacement[{a, b, c, d, e}, 10]

We can extract a range of elements using Take. Here we extract the first five elements of a list.

Take[{a, b, c, d, e, f, g}, 5]

Here we extract the last five elements.

Take[{a, b, c, d, e, f, g}, -5]

We can extract a sublist from the middle of the list, as well. Here we extract the third through sixth elements of a list.

Take[{a, b, c, d, e, f, g}, {3, 6}]

Drop works similarly, removing elements from a list. Here we drop the first three elements from a list.

Drop[{a, b, c, d, e, f, g}, 3]

Note that Take and Drop do not alter the value of a list. Here is a vector of values.

alpha
       = {𝒶, 𝒷, 𝒸, 𝒹, ℯ, 
        𝒻};

If we drop some elements from alpha, Mathematica returns the remaining elements.

Drop[alpha, {3, 5}]

However, the original list remains unchanged.

alpha

To alter the value of alpha, we must set the value of alpha to the desired subset.

alpha = Drop[alpha, {3, 5}];

alpha

Here is a function called sampleWithoutReplacement that takes a given number of samples from a list, without choosing the same element more than once.

sampleWithoutReplacement[L_List, n_Integer?Positive] := 
      Module[{ℒ = L, ℓ = Length[L], pos}, 
        Table[pos = Random[Integer, {1, ℓ}]; 
          ℒ = Drop[ℒ, {pos}]; 
          ℓ = ℓ - 1; L[[pos]], {n}]]

Here is a sample of fifteen data points from a list of integers from 1 to 50.

sampleWithoutReplacement[Range[50], 15]

Exercises: Extracting Parts of Lists

Using Part, extract the b from {{a, b, c}, {d, e, f}}.

The b is in the first row and second column of the matrix, so we use Part as follows.

Part[{{a, b, c}, {d, e, f}}, 1, 2]

We can use the double-bracket shorthand for Part as follows.

{{a, b, c}, {d, e, f}}[[1, 2]]

Switch the third and fourth elements of the list {a, b, c, d, e}.

We can use the part-selection notation [[{ }]] to switch the elements.

{a, b, c, d, e}[[{1, 2, 4, 3, 5}]]

Use Table and Part to define a function called everyOtherElement that extracts every other element of a list.

The key is to set up an iterator inside Table that extracts every second element from a list.

everyOtherElement[L_List] := Table[Part[L, n], {n, 1, Length[L], 2}]

everyOtherElement[{a, b, c, d, e, f, g, h, i, j, k, l}]

Define a function matrixTrace that computes the sum of all the elements along the diagonal of a matrix.

Here is a matrix.

testmat=(({{a, b, c}, {d, e, f}, {g, h, i}}));

The elements along the diagonal have equal coordinates in the horizontal and vertical directions.

Table[testmat[[n, n]], {n, 1, 3}]

To add the elements, we use Sum.

Sum[testmat[[n, n]], {n, 1, 3}]

Here is a general definition of matrixTrace.

matrixTrace[m_] := Sum[m[[n, n]], {n, Length[m]}]

Here is a square matrix.

primemat = Partition[Prime[Range[16]], 4]; MatrixForm[primemat]

The function works as expected.

matrixTrace[primemat]

A hand calculation shows that the answer is correct.

2 + 13 + 31 + 53

Here is a modification to allow only square matrices.

Clear[matrixTrace]

matrixTrace[m_] := 
      Sum[m[[n, n]], {n, Length[m]}] /; 
        MatrixQ[m] && First[Dimensions[m]] == Last[Dimensions[m]]

Characteristics of Elements

Often it is necessary to select elements from a list according to certain criteria, rather than according to position. The Select function returns the elements of a given list that meet a specified criterion.

For example, one way to extract all of the prime numbers less than 100 is first to generate a list of integers from 1 to 100, then select from that list all the elements that satisfy the predicate function PrimeQ.

Select[Range[100], PrimeQ]

Similarly we can choose from a list of ordered pairs all pairs that satisfy some specified criteria. Here is a list of points along a spiral.

pts = Table[{r Cos[r], r Sin[r]}, {r, 0, 25, 0.25}];

Here is a predicate function that acts on pairs of numbers that returns True if its first element times its second element is greater than zero.

quad1or3[{x_, y_}] := (x y > 0)

Here is a second predicate that returns True if its second element is greater than or equal to its first element.

ygex[{x_, y_}] := (y >= x)

(We want each of the following graphs to use the same option settings, and SetOptions allows us to change the default settings. We will change the settings back to their original values after this example.)

SetOptions[ListPlot, AspectRatio -> 1, 
      PlotRange -> {{-25, 25}, {-25, 25}}, 
      PlotStyle -> PointSize[0.02], DisplayFunction -> Identity, 
      Ticks -> None];

The first plot is of all the points in the list pts.

ptsplot = ListPlot[pts, PlotLabel -> "All points"];

The second plot is of the points in the first or third quadrant (those that satisfy the predicate function quad1or3).

q1or3pts = Select[pts, quad1or3];

q1or3plot = ListPlot[q1or3pts, PlotLabel -> "Quad. 1 or 3"];

The last plot is of the points for which the y value is greater than or equal to the x value.

ygexpts = Select[pts, ygex];

ygexplot = 
      ListPlot[ygexpts, PlotLabel -> "y ≥ x"];

Here are the three plots.

Show[
  GraphicsArray[{ptsplot, q1or3plot, ygexplot}]
    ];

Here we set the options to ListPlot back to the default settings.

SetOptions[ListPlot, AspectRatio -> (GoldenRatio^(-1)), 
      PlotRange -> Automatic, PlotStyle -> Automatic, 
      DisplayFunction :> $DisplayFunction, Ticks -> Automatic];

We can test whether an element is present in a list by using MemberQ and FreeQ.

MemberQ[{a, b, c, d, e}, c]

FreeQ[{a, (a^b), (a^((b^c)))}, c]

Exercises: Characteristics of Elements

Use Select to extract all the numbers from {a, 1, b, 2, c, 3}.

We use the built-in predicate function NumberQ.

Select[{a, 1, b, 2, c, 3}, NumberQ]

Create the following list of x, y pairs.

pts = 
      Flatten[Table[{x, y}, {x, -3, 3, 0.3}, {y, -3, 3, 0.3}], 1];

Select from it all points whose coordinates satisfy sin(x) < sin(y).

Here is a plot of all the points.

ListPlot[pts, AspectRatio -> 1, PlotStyle -> PointSize[0.02]];

Here is a predicate function that accepts an ordered pair and returns True if the sine of the horizontal coordinate is less than the sine of the vertical coordinate.

inRegionsQ[{x_, y_}] := Sin[x] < Sin[y]

We use Select to extract the points that satisfy inRegionsQ. The following is a shortened version of the output.

Short[Select[pts, inRegionsQ], 3]

Here is a plot of the points that satisfy inRegionsQ.

ListPlot[Select[pts, inRegionsQ], AspectRatio -> 1, 
      PlotStyle -> PointSize[0.02]];

Clear[pts, inRegionsQ]

Using Names and StringLength, find all Mathematica function names that contain more than twenty-five letters. Find all function names that contain fewer than three letters.

Here is a shortened list of all built-in Mathematica function names.

Short[Names["System`*"], 3]

Here is a predicate function that returns True if its argument is a function whose name contains more than twenty-five letters.

longNameQ[s_String] := StringLength[s] > 25

We use Select to extract the desired function names.

Select[Names["System`*"], longNameQ]

The following function returns True if its argument contains fewer than three letters.

shortNameQ[s_String] := StringLength[s] < 3

Here are all the names with one or two letters.

Select[Names["System`*"], shortNameQ]

List Structure

Characteristics of Structure

There are many functions to determine characteristics of a list, such as its length, depth, or dimensions; further, we can test if a list is a vector or matrix.

The length of a vector is the number of elements it contains.

vec = {a, b, c};

Length[vec]

Because the vector vec is one-dimensional, Dimensions[vec] returns a list containing the size of vec in its one dimension.

Dimensions[vec]

The length of a matrix is the number of rows it contains; that is, the number of sublists it contains.

Length[rect]

The dimensions of a matrix are its number of rows and columns.

Dimensions[rect]

In general, the length of an n-dimensional list is the number of n–1-dimensional sublists it contains. The list tensor contains two three-dimensional sublists, so its length is 2.

Length[tensor]

Dimensions returns the length in each “direction” of a nested list.

Dimensions[tensor]

If we have an irregularly shaped array of data, we can find the length of each sublist by mapping the function Length over the array.

irreg = {{1, 2, 3, 4}, {2, 4, 6}, {3, 6, 9, 12, 15, 18}, {4, 8}}; 
    TableForm[irreg]

Map[Length, irreg]

Here is a function that pads (with zeros) an irregularly shaped array, so that each sublist has the same length.

padArray[𝒜_List] := 
      Module[{𝓂𝒶𝓍 = 
            Max[Map[Length, 𝒜]]}, 
        Map[Take[Join[#, Table[0, {𝓂𝒶𝓍}]], 
              𝓂𝒶𝓍]&, 𝒜]]

padArray[irreg] // TableForm

Changing List Structure

Flatten is used to remove all but the outermost list brackets in a nested list of lists.

Flatten[hurr]

Here is a nested list called prods.

prods = Table[m×n, {m, 2, 5}, {n, 1, 5}]

Flatten removes all of the inner list brackets from prods, returning a one-dimensional list of values.

Flatten[prods]

The flattened list can now be sorted, have duplicates removed, and so forth.

Sort[Flatten[prods]]

Union[Flatten[prods]]

Lists can also be partitioned into sublists with a given length. Below a list of integers from 1 through 12 is partitioned into sublists of length 2, 3, and 4.

Partition[Range[12], 2]

Partition[Range[12], 3]

Partition[Range[12], 4]

Notice that all sublists created by Partition must be of the same length: if we partition a length-12 list into sublists of length 5, the two final elements are dropped.

Partition[Range[12], 5]

It is also possible to partition a list into “overlapping” sublists by giving the offset between subsequent lists as a third argument to Partition.

Partition[{a, b, c, d, e}, 3, 1]

Partition[{a, b, c, d, e}, 3, 2]

Often we have lists of data to be manipulated. Manipulating vectors (one-dimensional lists) is straightforward, but to manipulate multidimensional lists may require manipulation of the list’s structure.

For example, say we have a list of ordered pairs called mydata.

mydata = Table[{n, Sin[n + 0.1 Random[]]}, {n, 0, 9, 0.3}]

Transposing mydata results in a list containing two lists, the first being a list of x-coordinates, and the second being a list of y-coordinates.

Transpose[mydata]

With mydata in its transposed form, we can select the first sublist (the list of all x-coordinates) and manipulate it separately, and so forth. For example, here are the original points.

ListPlot[mydata, PlotJoined -> True];

Here we give the name tdata to the transpose of mydata.

tdata = Transpose[mydata]

Extracting the first part of the transpose of mydata returns a list containing all the x-values in mydata.

tdata[[1]]

Once the x-values have been extracted, we can manipulate them. Here we take the natural logarithm of 1 plus each of the elements.

Log[1 + %]

To get the transformed data back into the form of ordered pairs, we create a list containing a list of the transformed x-values, then a list of the original y-values.

transdata = {%, tdata[[2]]}

Transposing the previous list returns a list of ordered pairs.

final = Transpose[transdata]

Here is a graph of the transformed data.

ListPlot[final, PlotJoined -> True];

Notice that one effect of Transpose is to reverse the dimensions of a matrix.

Dimensions[mydata]

Dimensions[Transpose[mydata]]

Here are the dimensions of our array of hurricane data.

Dimensions[hurr]

Transposing the array returns a list containing a list of all the hurricane names, a list of all the wind speeds, and a list of all the pressures.

transhurr = Transpose[hurr]

Given the transformed list, we can obtain ordered pairs of wind speeds and pressures by first extracting the respective lists.

transhurr[[{2, 3}]]

Next we transpose the previous list to obtain the desired pairs.

hurrpairs = Transpose[transhurr[[{2, 3}]]]

The data are now in a form to be plotted using ListPlot.

ListPlot[hurrpairs, 
      FrameLabel -> {"Wind Speed", "Pressure"}, 
      PlotLabel -> 
        StyleForm["1995 Tropical Storm Data", 
          FontFamily -> "Helvetica", FontSize -> 10], 
      PlotStyle -> PointSize[0.02], Frame -> True, 
      RotateLabel -> True, TextStyle -> {FontSize -> 9}];

Exercises: List Structure

Define a function percentageOfPrimes that returns the percentage of prime numbers in a list.

To determine the percentage of primes in a list, we divide the number of primes by the length of the list. One way to determine the number of primes in a list is to use Length and Select.

Length[Select[Range[100], PrimeQ]]

A more straightforward method is to use Count and the conditional pattern _?PrimeQ.

Count[Range[100], _?PrimeQ]

Here is a definition of percentageOfPrimes.

percentageOfPrimes[𝒫_List] := 
      Length[Select[𝒫, PrimeQ]]/Length[𝒫]

Here is the percentage of primes in a list of integers from 1 to 100.

percentageOfPrimes[Range[100]]

Here is the percentage of primes in a list of odd numbers from 1 to 100.

percentageOfPrimes[Range[1, 100, 2]]

Generate the following list.

pts = Table[{r Cos[r], r Sin[r]}, {r, 0.0, 25.0, 0.25}];

Find the maximum x-value. Find the minimum y-value.

To extract the list of x-values, we take the first part of the transpose of pts. Once the list of x-values has been extracted, we use Max to determine the maximum.

Max[Transpose[pts][[1]]]

We follow a similar process to determine the minimum of the y-values.

Min[Transpose[pts][[2]]]

Sorting and Grouping Data

Sort sorts a list of numbers into increasing order.

Sort[{1, 3, 5, 7, 9, 2, 4, 6, 8}]

We can specify a function that determines if two adjacent elements are sorted. Here we tell Mathematica to sort a list of numbers in decreasing order, by specifying that two elements are sorted if they satisfy the function Greater (that is, the first element is greater than the second element).

Sort[{1, 3, 5, 7, 9, 2, 4, 6, 8}, Greater]

When working with multidimensional data, it is often necessary to explicitly define a sorting function. Here we read in the list of hurricane data.

Get["hurr.dat"]

To sort the list in ascending order of wind speed, we first define a sorting function as follows. The function returns True if the second element (the wind speed) of one record is less than the second element of another.

hurrSort[rec1_, rec2_] := rec1[[2]] < rec2[[2]]

Sort[hurr, hurrSort] // TableForm

Here is a pure sorting function equivalent to hurrSort.

Sort[hurr, #1[[2]] < #2[[2]]&] // TableForm

Split partitions a list into sublists that have identical elements. Here is a list of results from tossing a coin twenty times.

tosses = Table[If[Random[] < 0.5, "H", "T"], {20}]

Split divides the list into lists of runs.

runs = Split[tosses]

By mapping the function Length over the results of Split, we can compute the length of each run of heads or tails.

Map[Length, runs]

As an extension of the previous example, we generate a list containing pairs of the form (result, length of run).

Map[{First[#], Length[#]}&, runs]

This is called run length encoding.

Exercises: Sorting and Grouping Data

Using RotateRight, compute the first differences of a list of data.

Here is a list of data generated from a quadratic curve.

quaddata = Table[(n^2) - 5 n + 1, {n, 0, 10}]

We compute the first difference by rotating each element of the list of data one position to the right, and subtracting the rotated list from the original list. Here is the list of first differences.

quaddata - RotateRight[quaddata]

The first element of the subtracted list is meaningless (the first element of the original list minus the last element of the original list), so we use Drop to remove it.

Drop[%, 1]

Extract the five largest elements from a list of random numbers.

Here is a list of random numbers.

randomdata = Table[Random[], {15}]

Sort arranges the list into ascending order.

Sort[randomdata]

To extract the five largest elements in the original list, we extract the last five elements from the sorted list.

Take[%, -5]

Define a function called runEncode that takes a list of discrete random data and converts it into lists of the form (element, length of run). Define a function called runDecode that reconstructs a list resulting from runEncode.

Here is some coin-tossing data taken from the discussion.

tosses = Table[If[Random[] < 0.5, "H", "T"], {20}]

runEncode[ℒ_List] := 
      Map[{First[#], Length[#]}&, Split[ℒ]]

td = runEncode[tosses]

This definition of runDecode reverses the encoding of runEncode.

runDecode[ℒ_List] := 
      Flatten[Map[Table[First[#], {Last[#]}]&, ℒ]]

Here are the original data.

runDecode[td]

% == tosses

Graphing Lists of Data

Built-In Functions

There are many built-in functions for plotting lists of data. ListPlot accepts a list of ordered pairs, and generates a two-dimensional scatter plot of the points.

The file circle.dat defines a variable circledata to be a list of ordered pairs that fall on a circle.

<< circle.dat;

ListPlot[circledata];

Graphing functions take many options, which control how a graph is rendered. Here is the same graph with several options changed.

ListPlot[circledata, PlotJoined -> True, 
      GridLines -> Automatic, AspectRatio -> 1, 
      PlotLabel -> "circledata", AxesStyle -> Thickness[0.02]];

Options[graphing function] returns a list of options that a particular graphing function accepts, along with their default values.

ListPlot also accepts a flat list of data, and generates a scatter plot in which each data point is assumed to have horizontal coordinate 1, 2, 3, …. Here we create a scatter plot of the second coordinate of each ordered pair in circledata.

ListPlot[Last[Transpose[circledata]], PlotStyle -> PointSize[0.03]];

Clear[circledata]

The function ListPlot3D accepts a rectangular grid of height values, and generates a surface plot of the data. Here is a rectangular array of numbers.

newdata = Table[Binomial[n, m], {n, 1, 10}, {m, 1, 5}]; 
    TableForm[newdata]

Here is a surface plot of the array.

ListPlot3D[newdata];

Three-dimensional plotting commands accept options, as well.

ListPlot3D[newdata, ColorFunction -> Hue, 
      BoxRatios -> {1, 2, 0.5}, 
      ViewPoint -> {-2.5,  -1.5,  1.5}, 
      AxesLabel -> {n, m, Binomial[n, m]}, 
      PlotLabel -> "binomial data", Mesh -> False, 
      FormatType -> TraditionalForm];

Clear[newdata]

Similarly, we can create density plots and contour plots from lists of data, using ListDensityPlot and ListContourPlot. ListDensityPlot is useful for plotting digitized image data.

ListDensityPlot[
      ReadList["icondata.dat", Number, RecordLists -> True]];

Note that the built-in plotting commands expect a list of real numbers, and will print an error message if we try to plot a list containing any complex numbers, strings, or undefined symbols. For example, the package Miscellaneous`ChemicalElements` contains data related to the chemical elements.

Needs["Miscellaneous`ChemicalElements`"]

Here we create a scatter plot of the atomic weight divided by the atomic number of each element. (Mathematica warns that several of the elements do not have a stable isotope, which we can ignore for this example.)

elementdata = AtomicWeight[Elements]/AtomicNumber[Elements];

Here are the first ten data points of elementdata.

Take[elementdata, 10]

Here is a scatter plot of the ratios.

ListPlot[elementdata];

However, the melting point of an element is returned with the unit Kelvin, making the data unusable directly with ListPlot.

MeltingPoint[Hydrogen]

To remove the units from each measurement, we use algebraic division to divide out the unit symbols. (Some warning messages have been deleted.) Mathematica warns that the melting point of some elements is unknown.

elementdata = MeltingPoint[Elements]/Kelvin;

In the unknown cases, Mathematica returns the symbol Unknown, instead of a numeric value.

Last[elementdata]

When we use ListPlot to graph a list of points that contain the symbol Unknown, Mathematica prints a warning message that we are trying to plot nonnumeric data, and generates a graph of the remaining points.

ListPlot[elementdata, PlotJoined -> True];

Clear[elementdata]

Add-On Packages

The package Graphics`Graphics` defines several functions useful for plotting data. Here we load the package.

Needs["Graphics`Graphics`"]

The functions LogListPlot, LogLinearListPlot, and LogLogListPlot create various logarithmic plots of data.

wts = AtomicWeight[Elements];

LogListPlot[wts];

BarChart, GeneralizedBarChart, StackedBarChart, and PercentileBarChart create various bar charts and histograms of data, and PieChart creates pie charts.

Show[
      GraphicsArray[{
          BarChart[
            Transpose[{Take[wts, 5], Abbreviation[Take[Elements, 5]]}], 
            DisplayFunction -> Identity], 
          PieChart[Take[wts, 5], 
            PieLabels -> Abbreviation[Take[Elements, 5]], 
            DisplayFunction -> Identity]}]];

The various options to BarChart and PieChart provide many ways to customize the appearance of data graphs.

Clear[wts]

ErrorListPlot allows us to plot data along with error bars. Here are some points along a line with some noise added, each with an error bar 0.5 units long in each direction.

linedata = 
      Table[{x, x + Random[Real, {-0.5, 0.5}], 0.5}, {x, 0.0, 10.0, 
          0.5}];

ErrorListPlot[linedata, 
      Epilog -> {GrayLevel[0.5], Line[{{0, 0}, {10, 10}}]}];

Clear[linedata]

The package Graphics`MultipleListPlot` allows us to plot multiple lists of data on the same set of axes.

Needs["Graphics`MultipleListPlot`"]

Here is a scatter plot two lists of data.

MultipleListPlot[Log[MeltingPoint[Take[Elements, 10]]/Kelvin], 
      Log[BoilingPoint[Take[Elements, 10]]/Kelvin], PlotJoined -> True];

By changing the options to MultipleListPlot, we can change the symbol used to plot each data set, include error bars or ellipses, and much more.

Strings

In many cases it is convenient to keep data in the form of a string, and to manipulate the string directly. Mathematica contains many string manipulation functions that mirror list manipulation commands.

Reverse[{b, a, c, k, w, a, r, d, s}]

StringReverse["backwards"]

Join[{a, b, c}, {x, y, z}]

StringJoin["abc", "xyz"]

The file register.dat contains the following data.

!! register.dat

The first eight columns of each line contain the name of a state, the ninth through twelfth columns contain the number of automobiles registered (in thousands), the thirteen through sixteenth columns contain the number of thousands of trucks and buses registered, the seventeenth through twentieth columns contain the number of thousands of motorcycles registered, and the remaining columns contain the total for each state. (The first line of data is used to number the columns.)

Here we read in each line as a string, discarding the first line of data.

register = Rest[ReadList["register.dat", String]]

The following function dissects each string into its respective list of data.

dissect[𝒮_String] := {StringTake[𝒮, 8], 
        StringTake[𝒮, {9, 12}], 
        StringTake[𝒮, {13, 16}], 
        StringTake[𝒮, {17, 20}], 
        StringTake[𝒮, {21, 24}]}

Map[dissect, register] // TableForm

As defined above, dissect leaves the numbers in each measurement in the form of a string. The following addition to the definition of dissect changes the numerical strings into numbers.

dissect[𝒮_String] := 
      Flatten[{StringTake[𝒮, 8], 
          ToExpression[{StringTake[𝒮, {9, 12}], 
              StringTake[𝒮, {13, 16}], 
              StringTake[𝒮, {17, 20}], 
              StringTake[𝒮, {21, 24}]}]}]

Map[dissect, register] // TableForm

Just as parts of lists and expressions can be replaced using ReplaceAll (/.), parts of strings can be replaced using StringReplace. In the following string we replace every occurrence of the word “Sin” with “sin” and replace every square bracket with a matching parenthesis.

StringReplace[
      "Sin[Sin[x]]", {"Sin" -> "sin", 
        "[" -> "(", "]" -> ")"}]

Exercises: Strings

How many characters are in the word Mathematica?

StringLength returns the number of characters in a string.

StringLength["Mathematica"]

Using Characters and Permutations, find all anagrams (rearrangements of letters) of the word math.

Here are the characters in the string “math”.

Characters["math"]

Permutations returns all possible permutations of the elements in a list.

anagrams = Permutations[Characters["math"]]

There are 24 (4!) permutations.

Length[%]

To combine each list of characters into a string, we use StringJoin.

Table[StringJoin[anagrams[[i]]], {i, 24}]

Alternatively, we can use the high-level programming function Map to map StringJoin over each list of characters.

Map[StringJoin, anagrams]

3. Curve Fitting and Interpolation

Linear Regression

Fit

The Fit command creates a best-fitting function, using a given list of basis functions, for a set of points. Here is a list of ordered pairs.

data  =  {{1, 4}, {2, 3}, {4, 3}, {7, 6}, {10, 11}};

Here is a plot of the data.

dataplot = 
      ListPlot[data, PlotStyle -> {PointSize[0.03], GrayLevel[0.5]}];

We can fit a quadratic function of the form a + b x + c (x^2) (that is, find values for the parameters a, b, and c that produce the best-fit curve) to the data set by giving Fit a list of basis functions to use ({1, x, x^2} in this case), and the variable in which to create the polynomial.

fitcurve = Fit[data, {1, x, (x^2)}, x]

Here is a plot of the best-fit curve.

curveplot = Plot[fitcurve, {x, 1, 10}, PlotStyle -> Hue[0.05]];

Here the data plot and the best-fit plot are combined.

Show[{dataplot, curveplot}];

To find a value of the fit function at a particular point, we use a replacement rule. Here is the value of the fit function at x = 6.

fitcurve /. x -> 6

We can fit the same data using different sets of basis functions. Here is a best-fitting function made up of a linear combination of a constant, an exponential, and a logarithm function.

fit2[x_] = Fit[data, {1, Exp[x], Log[x]}, x]

Here we use the Epilog option of Plot to have Mathematica plot the data points on the same axes as the plot.

Plot[fit2[x], {x, 1, 10}, PlotRange -> {2, 12}, 
      PlotStyle -> Hue[0.6], 
      Epilog -> {PointSize[0.03], GrayLevel[0.5], Map[Point, data]}];

Clear[data, dataplot, fitcurve, curveplot, fit2]

There are two cases when the Fit command is inappropriate. The first occurs when the function which we want to fit to is not linear in the unknown parameters. In such cases, it is necessary to use the package Statistics`NonlinearFit`, described in the next section. The second case occurs when the data set contains outlying points which a least-squares fit will weight too heavily. In such cases, a more robust fitting process should be used. A more rigorous technique is defined in the add-on Experimental Data Analyst.

Note that Fit is not limited to two-dimensional data. Here is a list of ordered triples.

wavydata = 
      Flatten[Table[{i, j, 
            2  (i^2) - Cos[j] + Random[Real, {-0.2, 0.2}]}, {i, -3.0, 
            3.0, 0.2}, {j, -6.0, 6.0, 0.2}], 1];

Here is the first ordered triple.

First[wavydata]

The function ScatterPlot3D, defined in the package Graphics`Graphics3D`, displays a list of points in three dimensions.

Needs["Graphics`Graphics3D`"]

wavydataplot = 
      ScatterPlot3D[wavydata, BoxRatios -> {1, 1, 0.4}];

Here we give Fit a list of basis functions in two variables. The result is a linear combination of the basis functions.

wavefit = Fit[wavydata, {1, (x^2), Cos[y]}, {x, y}]

Again, we can compute a value of the fit function for particular values of its parameters by using replacement rules.

wavefit /. {x -> 2, y -> -3}

Here is a plot of the best-fit surface.

wave = Plot3D[wavefit, {x, -3, 3}, {y, -6, 6}];

Here are the data and best-fit surface combined.

Show[{wave, wavydataplot}];

Clear[wavydata, wavydataplot, wavefit, wave]

Regress

The standard package Statistics`LinearRegression` defines a function called Regress, which allows us to extract more information during the fitting process. Here we load the package.

Needs["Statistics`LinearRegression`"]

Regress[data, funs, vars] finds a least-squares fit to a list of data as a linear combination of the functions funs, in the variables vars. The data can have the form {{(Subscript[x, 1]), (Subscript[y, 1]), …, (Subscript[f, 1])}, {(Subscript[x, 2]), (Subscript[y, 2]), …, (Subscript[f, 2])}, …}, where the number of coordinates x, y, … is equal to the number of variables in the list vars. The data can also be of the form {(Subscript[f, 1]), (Subscript[f, 2]), …}, with a single coordinate assumed to take values 1, 2, …. The argument funs can be any list of functions that depend only on the objects vars.

Here are some random data points.

cosdata  =  
      Table[{x, Cos[x] + Random[Real, {-0.1, 0.1}]}, {x, -2.0, 2.0, 
          0.1}];

ListPlot[cosdata, PlotStyle -> {PointSize[0.015], GrayLevel[0.5]}];

Fit returns only the best-fit function.

Fit[cosdata, {1, x, (x^2), (x^3)}, x]

Regress, defined in the package Statistics`LinearRegression`, by default returns various test statistics, including the coefficient of determination ((R^2) value), adjusted coefficient of variation (\(TraditionalForm\`R\&_\^2\) value), and an ANOVA table.

Regress[cosdata, {1, x, (x^2), (x^3)}, x]

In the above parameter table, we see that the p-values for the coefficients of x and (x^3) are large, which suggests that the “true” values of the parameters are zero. Here is a report based on a model containing only a constant and quadratic term.

Regress[cosdata, {1, (x^2)}, x]

We can compute different statistics for a least-squares fit by changing the value of the option RegressionReport. A list of all possible values for RegressionReport is returned when we enter the following.

RegressionReportValues[Regress]

Here we create a report that returns only the best-fit function, the predicted values, and the residual values.

report = 
      Regress[cosdata, {1, (x^2)}, x, 
        RegressionReport -> {BestFit, PredictedResponse , FitResiduals}];

Here we extract the best-fit function, using a replacement rule.

BestFit /. report

Here are abbreviated versions of the lists of predicted responses and fit residuals.

Short[PredictedResponse /. report, 2]

Short[FitResiduals /. report, 3]

Here is a plot of the residuals.

residualplot = 
      ListPlot[FitResiduals /. report, FrameTicks -> None, 
        Frame -> True, Ticks -> None, PlotStyle -> Hue[0], 
        TextStyle -> {FontSize -> 9}, 
        PlotLabel -> "residuals"];

Here we combine several results on the same graph.

Plot[Evaluate[BestFit /. report], {x, -2, 2}, 
      Epilog -> {Rectangle[{1, 0.6}, {2.1, 1}, residualplot], 
          PointSize[0.01], GrayLevel[0.5], Map[Point, cosdata]}, 
      PlotLabel -> "Quadratic fit"];

Clear[cosdata, report, residualplot]

Exercises: Linear Regression

Read in the contents of the file primes.dat, and fit a line and a quadratic function to the data.

Here we read in the data file.

primedata = ReadList["primes.dat", Number]

Here is a plot of the data.

primeplot = ListPlot[primedata, PlotStyle -> PointSize[0.02]];

Here is a line fit to the data, using the basis functions 1 and x.

linfit = Fit[primedata, {1, x}, x]

Here is a quadratic curve fit to the data, using the basis functions 1, x, and (x^2).

quadfit = Fit[primedata, {1, x, (x^2)}, x]

Here are the two curves.

Plot[{linfit, quadfit}, {x, 0, 17}, 
      PlotStyle -> {GrayLevel[0.5], GrayLevel[0.25]}];

Here we combine the curves with the original data.

Show[{%, primeplot}];

Clear[primedata, primeplot, linfit, quadfit]

Nonlinear Fitting

When we need to fit data to a model that is not a linear combination of basis functions, the built-in function Fit cannot be used. For example, here are some data that are generated by the function 0.75 sin(1.5 x), with some noise added.

sindata = 
      Table[{x, 0.75 Sin[1.5 x] + Random[Real, {-0.1, 0.1}]}, {x, 
          -3.0, 3.0, 0.1}];

Here is a plot of the data.

ListPlot[sindata];

The model for the data, a sin(b x), is not a linear combination of basis functions, so it is necessary to use the function NonlinearFit from the package Statistics`NonlinearFit`. Here we load the package.

Needs["Statistics`NonlinearFit`"]

NonlinearFit[data, model, vars, params] uses iterative methods to fit the data to the model containing the given parameters. Here is the NonlinearFit command used for our data.

sinfit = NonlinearFit[sindata, a Sin[b x], x, {a, b}]

Here is a plot of the model generated by Mathematica, shown with the original data.

Plot[sinfit, {x, -3, 3}, PlotRange -> All, 
      Epilog -> {PointSize[0.02], RGBColor[0, 0, 1], 
          Map[Point, sindata]}];

Clear[sindata, sinfit]

Spline Fitting

The package NumericalMath`SplineFit` allows us to fit a spline to data for which the fit is not one-to-one; that is, we wish to approximate data with a general curve in the plane, rather than a one-to-one function. Here we load the package.

Needs["NumericalMath`SplineFit`"]

Here we read into Mathematica a list of corners that make up the corners of a map of the United States. (Similar data for all countries is found in the add-on package Miscellaneous`WorldData`.)

Needs["Miscellaneous`WorldData`"]

Here we extract the data for the continental United States, numericize the data, and then reverse each (latitude, longitude) pair in order to get the proper (x, y) shape for graphing.

USdata = Map[Reverse, N[First[WorldData["USA"]]]];

The data are in the form of (x, y) pairs.

First[USdata]

Here is a plot of our data points.

ListPlot[USdata, AspectRatio -> Automatic, 
      PlotJoined -> True, PlotStyle -> RGBColor[0, 1, 0]];

Possible types of splines we can use are cubic, Bézier, and composite Bézier. Here we generate a cubic spline based on the data.

spline = SplineFit[USdata, Cubic]

The result is a SplineFunction object that displays the type of spline and the range of the single parametric argument. (The domain of the parameter is from zero to one less than the number of data points.) To compute the (x, y) value corresponding to a value of the parameter, we give the SplineFunction object a value inside square brackets.

spline[10]

Here is a plot of the spline function, shown together with the data points. We see that a cubic spline interpolates all the data points.

ParametricPlot[spline[t], {t, 0, 83}, Compiled -> False, 
      AspectRatio -> Automatic, 
      Prolog -> {RGBColor[0, 1, 0], PointSize[0.015], 
          Map[Point, USdata], GrayLevel[0]}];

Here is a Bézier spline based on the data. A Bézier spline interpolates only the endpoints of a data set, while the other points control the shape of the spline.

ParametricPlot[
      Evaluate[SplineFit[USdata, Bezier]][t], {t, 0, 83}, 
      Compiled -> False, AspectRatio -> Automatic, 
      Prolog -> {RGBColor[0, 1, 0], PointSize[0.015], 
          Map[Point, USdata], GrayLevel[0]}];

A composite Bézier spline interpolates every other data point.

ParametricPlot[
      Evaluate[SplineFit[USdata, CompositeBezier]][t], {t, 0, 83}, 
      Compiled -> False, AspectRatio -> Automatic, 
      Prolog -> {RGBColor[0, 1, 0], PointSize[0.015], 
          Map[Point, USdata], GrayLevel[0]}];

Interpolation

Basic Interpolation

In the file lifedata.dat we have data taken from a life table, stored in a variable called lifedata. Here we execute the file.

Get["lifedata.dat"]

We have measurements of the number of survivors in a group of 100,000 persons, taken ten years apart. Using interpolation, we can define a function which will take this data set and return a continuous function that estimates measurements at any time in the domain of the data. The function Interpolation returns an InterpolatingFunction object, which is a continuous function over the domain of the data. Here we add [x] to the right-hand side of the assignment to create a function which evaluates the InterpolatingFunction at the point x.

life[x_] = Interpolation[lifedata][x]

This function evaluates to a known value when evaluated at one of the data points.

life[50]

Elsewhere, the data are interpolated to provide an estimate.

life[55]

Next we plot the data points and the interpolation function simultaneously.

Plot[life[x], {x, 0, 60}, 
      Epilog -> {PointSize[0.03], Map[Point, lifedata]}];

Interpolation has one option, InterpolationOrder, which specifies the order of the polynomial used to join each point to its nearest neighbors. The default value is 3, meaning that the curve between any two points is obtained by fitting a cubic equation to these points and their nearest neighbors. Setting InterpolationOrder to one we obtain piecewise linear interpolation.

lininterp[x_] = 
      Interpolation[lifedata, InterpolationOrder -> 1][x]

Here is a plot of our piecewise linear interpolating function.

Plot[lininterp[x], {x, 0, 60}, 
      Epilog -> {PointSize[0.03], Map[Point, lifedata]}];

Here we create interpolating functions of different orders.

ordertable[x_] = 
      Table[Interpolation[lifedata, InterpolationOrder -> a][x], {a, 
          1, 4}];

Here is a plot of the interpolating functions over part of the domain of the data.

Plot[Evaluate[ordertable[x]], {x, 0, 12}, 
      Epilog -> {PointSize[0.03], Map[Point, lifedata]}, 
      PlotStyle -> Table[Hue[i/4], {i, 4}]];

Mathematica allows us to differentiate interpolating functions. Because interpolating functions are piecewise functions, the derivative may not be a smooth function. Here is the derivative of our piecewise linear interpolating function.

Plot[Evaluate[lininterp]'[t], {t, 0, 60}];

Mathematica allows us to extrapolate values that fall outside the domain of the data set, printing a warning message that extrapolation is being used.

life[65]

List Interpolation

To construct an approximate function that represents an arbitrary array of data, we use ListInterpolation. Here we reevaluate the file that sets eye to be an array of gray-level data taken from a photograph.

Get["eye.dat"];

Here we extract a submatrix of data, calling it eyedata.

eyedata = Map[Take[#, {40, 80}]&, Take[eye, {20, 40}]];

Here is a contour plot of the data set.

ListContourPlot[eyedata, AspectRatio -> Automatic, 
      ContourLines -> False];

Using ListInterpolation, we generate a continuous function that interpolates the data.

interpeye = ListInterpolation[eyedata]

Here is a contour plot of the interpolating function. Because the function is continuous, we can generate a smoother graph than is possible using the original data points.

ContourPlot[interpeye[x, y], {y, 1, 41}, {x, 1, 21}, 
      AspectRatio -> Automatic, PlotPoints -> {20, 40}, 
      ContourLines -> False];

Function Interpolation

FunctionInterpolation generates an interpolating function that interpolates values of a function over a given domain. Here we create an interpolating function that computes approximate values of a function for values of its parameter between -3 and 3.

approxNormalPDF = 
      FunctionInterpolation[
        (1/(Sqrt[2  π])) Exp[((-(t^2))/2)], {t, -3, 3}]

Here is a particular value of the approximate function.

approxNormalPDF[1.5]

Here is a plot of the function.

Plot[approxNormalPDF[x], {x, -3, 3}];

FunctionInterpolation also combines combinations of interpolating functions into a single InterpolatingFunction object. Here we create a single two-dimensional interpolating function from a product of interpolating functions.

approxNormProduct = 
      FunctionInterpolation[
        approxNormalPDF[x]× approxNormalPDF[y], {x, -3, 3}, {y, 
          -3, 3}]

Here is a three-dimensional plot of the new interpolating function over a square region.

Plot3D[approxNormProduct[s, t], {s, -1, 2}, {t, -2, 1}, 
      ColorFunction -> GrayLevel];

Exercises: Interpolation

Read in the contents of the file primes.dat and create a function that interpolates the data. Compute the values of the function corresponding to parameter values of 10 and 10.25.

primedata = ReadList["primes.dat", Number]

primedataInt = Interpolation[primedata]

The interpolating function exactly equals the tenth data point when given the argument 10.

primedataInt[10]

primedata[[10]]

The function returns an interpolated value at the argument 10.25.

primedataInt[10.25]

The file icondata.dat contains a 32 × 32 array of gray levels in a picture of a Mathematica notebook icon. Read in the file and create a two-dimensional interpolating function from the data.

We read in the data using ReadList.

icon = 
      ReadList["icondata.dat", Number, RecordLists -> True];

Here is the fifteenth row of gray levels in the array.

icon[[15]]

Here is a density plot of the data.

ListDensityPlot[icon, Mesh -> False];

ListInterpolation creates a two-dimensional interpolating function from the data.

iconInterpolation = ListInterpolation[icon]

Here is a plot of the interpolating function that uses more sample points than there are data points.

DensityPlot[
      Evaluate[iconInterpolation[s, t]], {t, 1, 32}, {s, 1, 32}, 
      PlotPoints -> 65, Mesh -> False];

We can manipulate the interpolating function just as we manipluate any other function.

Plot3D[Evaluate[1 - iconInterpolation[s, t]], {t, 5, 27}, {s, 5, 27}, 
      PlotPoints -> 45, Mesh -> False, ColorFunction -> Hue, 
      BoxRatios -> {1, 1, 0.1}, Axes -> {Automatic, Automatic, False}];

4. Statistical Analysis

Data Manipulation

The package Statistics`DataManipulation` defines many functions for manipulating lists of data. Here we load the package.

Needs["Statistics`DataManipulation`"]

DropNonNumeric[data] drops rows of data that contain nonnumeric elements. Here we drop the nonnumeric header row from a list of data.

DropNonNumeric[(({{"atomic no.", "atomic wt."}, {1, 1.00794}, {3, 6.941}, {11, 22.98977}}))]//TableForm

Similarly, DropNonNumericColumn[data] drops columns of data that contain nonnumeric elements.

CumulativeSums[list] gives successive sums of the elements of list.

CumulativeSums[ 1/Range[10] ]

Ordering[list] gives the permutation that puts the elements of list in order. QuantileForm[list] sorts the elements in list, then gives a list of the elements, paired with their quantile position.

Ordering[ {2, 3, 54, 1, 3, 5}]

QuantileForm[{2, 3, 54, 1, 3, 5}]

TakeWhile[list, pred] takes elements from the beginning of list while pred returns True. LengthWhile[list, pred] gives the number of elements at the beginning of list for which pred returns True. Here we create a list of numbers of the form (2^((2^n))) + 1. (These numbers are called Fermat primes.)""

fermatPrimes = (2^((2^Range)[0, 6])) + 1

Here we take the elements from the beginning of the list which satisfy PrimeQ.

TakeWhile[ fermatPrimes,  PrimeQ ]

Here we see that there are five consecutive prime numbers at the beginning of fermatPrimes.

LengthWhile[ fermatPrimes,  PrimeQ ]

Column[data, n] extracts the (n^th) column in data. Column[data, {(Subscript[n, 1]), (Subscript[n, 2]), …}] gives a list of columns (Subscript[n, 1]), (Subscript[n, 2]), ….

Column[(({{a, 1}, {b, 2}, {c, 3}, {d, 4}})),2]

ColumnTake[data, spec] takes the specified columns in data, using the same specifications as Take. ColumnDrop[data, spec] drops the specified columns in data. Here is a list of data.

data = Table[(1/(i + j - 1)), {i, 3}, {j, 6}]; TableForm[data]

Here we take the first two columns of the list.

ColumnTake[data,  2] // TableForm

Here we drop the second and third columns of the list.

ColumnDrop[data, {2, 3}] // TableForm

Grouping Data

Some common operations for grouping data are also available in the package Statistics`DataManipulation`. (We enter Needs["Statistics`DataManipulation`"] if the package has not been loaded.)

To illustrate these functions we first create a list of random real numbers.

randdata = Table[Random[Real, {0, 10}], {1000}];

ListPlot[randdata];

A common operation that allows us to better understand the distribution of the data is to group or “bin” the data. BinCounts[data, {xmin, xmax, dx}] gives a list of the number of elements in data that lie in bins from xmin to xmax in steps of dx. BinCounts works analogously with multidimensional data.

Because the data are uniformly distributed, the bin counts are approximately equal.

counts = BinCounts[randdata, {0, 10, 1}]

BarChart provides a convenient way to view the results.

Needs["Graphics`Graphics`"]

BarChart[counts];

RangeCounts allows us to specify a detailed description of the intervals. We consider one large bin from zero to five, and then five others from five to ten.

RangeCounts[randdata, {0, 5, 6, 7, 8, 9, 10}]

BinLists and RangeLists return arrays containing the data partitioned into bins instead of just the sizes of the bins. In the following example we sort the data into evenly spaced bins of size 0.1, and look at the contents of the first bin.

BinLists[randdata, {0, 10, 0.1}]  //  First

CategoryCounts and CategoryLists are generalizations of BinCounts and BinLists that allow us to extend this analysis to arbitrary elements and not just numbers. CategoryCounts[{(Subscript[x, 1]), (Subscript[x, 2]), …}, {(Subscript[e, 1]), (Subscript[e, 2]), …}] gives a list containing the number of elements in the data {(Subscript[x, 1]), (Subscript[x, 2]), …} that match each of the (Subscript[e, i]).

These functions are appropriate when the data take continuous values. If the data take only a finite number of discrete values, and we wish to count the number of times each value is attained, then we use the Frequencies function. We illustrate its use by simulating one thousand rolls of a pair of six-sided dice.

rolldie := Random[Integer, {1, 6}]

diedata = Table[rolldie + rolldie, {1000}];

Frequencies[data] returns a list of pairs of the form {frequency, value}.

diefreqdata = Frequencies[diedata]

Again, BarChart allows us to visualize the results.

BarChart[diefreqdata, PlotRange -> All];

Clear[randdata, counts, diedata, diefreqdata]

Descriptive Statistics

Statistics`DescriptiveStatistics` and Statistics`MultiDescriptiveStatistics`define many useful functions for describing data. The functions are defined to work with the statistical distributions defined in other packages. Here is a listing of all the symbols defined in these packages.

Needs["Statistics`DescriptiveStatistics`"]

?Statistics`DescriptiveStatistics`*

Here is a list of random data.

sampledata = Table[Random[Real, {0, 10}], {50}];

sampleplot = 
      ListPlot[sampledata, 
        PlotStyle -> {PointSize[0.02], RGBColor[0, 0, 1]}];

Descriptive statistics can be divided into three categories: location statistics, dispersion statistics, and shape statistics. Location statistics provide basic information on the location of the data. These functions are generally straightforward and their usage can be learned from their usage message. Note that for each of these categories there is a report-generating function (LocationReport, DispersionReport, and ShapeReport) that returns a list of several of the most common statistics.

sampleμ = Mean[sampledata]

samplesd = StandardDeviation[sampledata]

GeometricMean[sampledata]

Median[sampledata]

Quantile[sampledata, (1/9)]

Here we add a band of color denoting one standard deviation on either side of the mean.

\(\(Show[sampleplot, 
	
      Prolog -> {Hue[0], 
          Rectangle[{0, sampleμ - samplesd}, {50, 
              sampleμ + samplesd}]}]; \)\)

TrimmedMean[list, f] gives the mean of the entries in list, with a fraction f of entries at each end dropped.

TrimmedMean[sampledata, (1/8)]

LocationReport[sampledata]

Dispersion statistics give information on how spread out the data are. Since our data are uniformly distributed random numbers, we expect a relatively high variance and standard deviation.

DispersionReport[sampledata]

Some of the available shape statistics are Skewness, which describes the amount of asymmetry in the distribution of the data, and Kurtosis, which measures the concentration of data around the mean. Here is the shape report for our random data.

ShapeReport[sampledata]

The function ZeroMean transforms a data set to have zero mean, and Standardize transforms a data set to have zero mean and variance one.

Data Smoothing

The package Statistics`DataSmoothing` defines the functions MovingAverage, MovingMedian, ExponentialSmoothing, and LinearFilter. Here we load the package.

Needs["Statistics`DataSmoothing`"]

A Random Walk

Here we generate a list of points using a random walk. At each point, the next point is generated by adding a random real number between -0.5 and 0.5. We will then use the tools in Statistics`DataSmoothing` to smooth out the random nature of the walk.

addrand[orig_] := orig + Random[Real, {-0.5, 0.5}]

randomwalk = NestList[addrand, 0, 1000];

Here is one way to create a list of the form {i, (Subscript[randomwalk, i])}.

modwalk = Table[{i, randomwalk[[i]]}, {i, Length[randomwalk]}];

An alternative way to create the same list is the following.

modwalk = Transpose[{Range[Length[randomwalk]], randomwalk}];

Here is a list of graphics directives and primitives that will be used for comparison in later graphs.

modpoints = {Hue[0], PointSize[0.005], Map[Point, modwalk]};

Here is a plot of the points along the random walk.

Show[Graphics[modpoints], Axes -> True];

In the following sections we will examine various techniques for smoothing the data set describing the random walk.

Moving Averages

MovingAverage[datalist, n] is an n-term moving average.

A moving average is a list of means of a series of subsets of a data set. For example, the (100^th) moving average of a list of 1000 data points is a list with the following form.

{mean[(Subscript[i, 1]),…,(Subscript[i, 100])],  mean[(Subscript[i, 2]),…, (Subscript[i, 101])], …, mean[(Subscript[i, 901]),…, (Subscript[i, 1000])]}

We can smooth this random walk using MovingAverage. Adding the semicolon after the following statement tells Mathematica to compute the MovingAverage and assign the resulting list to the variable, but to suppress printing the result.

averaged = MovingAverage[randomwalk, 100];

We may now plot the result, and compare it with the original series. However, the original data and our averaged list are now of two different sizes, and we must fix the origin. We fix the origin by creating a list of x-values for the averaged data such that each point lies midway through the horizontal interval over which it was averaged.

Here we create a list of the form {51, 52, 53, …, 950} and combine it with averaged, which has the form {i1, i2, i3, …, i900}, to create a list of pairs of the form {{51, i1}, {52, i2}, …, {950, i900}}.

modavg = Table[{i + 50, averaged[[i]]}, {i, Length[averaged]}];

ListPlot[modavg, PlotJoined -> True, 
      PlotStyle -> Thickness[0.002], Epilog -> modpoints];

Moving Median

MovingMedian[datalist, r] smoothes datalist using a moving median of span r. Each element in the result is the median of r elements from datalist.

movmed = MovingMedian[randomwalk, 100];

modmm  =  
      Table[{i + 50,  movmed[[i]]},  {i, Length[movmed]}];

ListPlot[modmm, PlotJoined -> True, 
      PlotStyle -> Thickness[0.002], Epilog -> modpoints];

Linear Filter

LinearFilter[datalist, {(Subscript[c, 0]), (Subscript[c, 1]), (Subscript[c, 2]), …}] passes the data in datalist through a linear filter with weights (Subscript[c, 0]), (Subscript[c, 1]), (Subscript[c, 2]), …. (A moving average is a special case of a linear filter.) Each element in the resulting list is given by Sum[(Subscript[c, j])(Subscript[datalist, t - j]), {j, 0, r - 1}], for t = r, r + 1, …, n, where n is equal to the length of datalist, and r is the number of weights.

linfilter  =  LinearFilter[ randomwalk,  {0.5, 0.3, 0.1, 0.1}];

ListPlot[linfilter, PlotJoined -> True, 
      PlotStyle -> Thickness[0.002], Epilog -> modpoints];

Exponential Smoothing

ExponentialSmoothing[datalist, a] smooths datalist using an exponentially weighted average with smoothing constant a, taking the first entry of datalist as the starting value. ExponentialSmoothing[datalist, a,(Subscript[y, 0])] smooths datalist, taking the starting value to be (Subscript[y, 0]). If (Subscript[y, t]) is the smoothed result at time t, then the result at time t + 1 is given by (Subscript[y, t])+a ((Subscript[datalist, t + 1])-(Subscript[y, t])), where a is the smoothing constant (0 < a < 1).

expsmooth = ExponentialSmoothing[randomwalk, 0.05];

ListPlot[expsmooth, PlotJoined -> True, 
      PlotStyle -> Thickness[0.002], Epilog -> modpoints];

Here we clear the data variables.

Clear[randomwalk, modwalk, modpoints, averaged, modavg, movmed, modmm, 
      linfilter, expsmooth]

Fourier Transforms

A frequent use of Fourier transforms is to filter noisy signals. The noise is generally of high frequency and is easily removed when we have a description of the signal in frequency space. First we generate a fairly noisy signal.

noisysound = 
      Table[Cos[x] + Sin[3 x] + Random[Real, {-0.2, 0.2}], {x, 0.0, 6.0, 
          0.01}];

ListPlot[noisysound, PlotJoined -> True];

Fourier[list] finds the discrete Fourier transform of a list of numbers. We can look at the profile of this function in frequency space by taking the Fourier transform and using ListPlot to view the absolute values of the complex numbers in this list.

ListPlot[Abs[Fourier[noisysound]], PlotJoined -> True, 
      PlotRange -> {0, 1}];

Chop[number, tolerance] replaces any number that is within tolerance of zero with an exact zero. A very simple filtering procedure is to chop any element in the Fourier transform whose absolute value is less than a certain threshold (in this example, 0.5).

filtereddata = Chop[Fourier[noisysound], 0.5];

Here is a plot of the smoothed frequency profile.

ListPlot[Abs[filtereddata], PlotJoined -> True, 
      PlotRange -> {0, 1}];

We perform an inverse Fourier transform on the filtered data, then graph it to see the underlying function within the original data.

ListPlot[Chop[InverseFourier[filtereddata]], 
      PlotJoined -> True];

Statistical Distributions

Available Distributions and Related Functions

Mathematica has standardized descriptions of most continuous and discrete distributions. For each distribution, Mathematica knows the probability density function, its cumulative distribution function, how to compute quantiles, the domain of the variables, mean, variance, standard deviation, skewness, kurtosis; the characteristic function; and how to create random samples with that distribution. We will need definitions from several packages. In practice, it is easier to load the package Statistics`Master`, but here we load each package individually.

Here we load the package containing continuous statistical distributions.

Needs["Statistics`ContinuousDistributions`"]

These are the distributions defined in the package.

Names["Statistics`ContinuousDistributions`*Distribution"]

Here we load the package containing discrete statistical distributions, and list the distributions defined in it.

\(Needs["Statistics`DiscreteDistributions`"]; 

    Names["Statistics`DiscreteDistributions`*Distribution"]\)

As each of these packages are loaded, Mathematica automatically loads related common packages defining utility functions, including CDF, which gives the cumulative distribution function (c.d.f.); PDF, which gives the probability density function (p.d.f.); CharacteristicFunction, which returns the characteristic function for the specified distribution; Domain, which gives the domain of the specified distribution; and QuantileQ, which returns True if the argument is a valid quantile specification (a number between 0 and 1).

Here is the c.d.f. of the standard normal distribution.

CDF[NormalDistribution[0, 1], x]

Here is the p.d.f. of the same distribution.

PDF[ NormalDistribution[0, 1],  x]

Here is the characteristic function of a chi-square distribution with parameter 3.

CharacteristicFunction[ChiSquareDistribution[3], t]

Parameters for a distribution do not have to be numeric. Here is the variance of a binomial distribution with parameters n and p.

Variance[BinomialDistribution[n,  p]]

Using Distributions

NormalDistribution[μ, σ] represents the normal (Gaussian) distribution with mean μ and standard deviation σ. Here are two functions corresponding to a particular normal distribution.

CharacteristicFunction[NormalDistribution[0.5, 0.2], x]

CDF[NormalDistribution[0.5, 0.2], x]

Here we define the function normalpdf to be the p.d.f. of a normal distribution with mean 0.5 and standard deviation 0.2.

normalpdf[x_] = PDF[NormalDistribution[0.5, 0.2], x]

Here is a plot of normalpdf.

Plot[normalpdf[x], {x, -0.2, 1.2}];

We can use the statistical distributions defined in the statistics packages as models for fitting data. Here is the p.d.f. of the standard normal distribution.

dist = PDF[NormalDistribution[0, 1], x]

Here is a table of points that follow the standard normal distribution, with some noise added.

normdata = 
      Table[With[{i = Random[Real, {-3, 3}]}, {i, 
            Random[Real, {-0.05, 0.05}] + dist /. x -> i}], {100}];

Here is a plot of the data.

ListPlot[normdata];

We use the p.d.f. of a general normal distribution (with parameters μ and σ) as a model for use in the NonlinearFit function defined in the package Statistics`NonlinearFit`.

model = PDF[NormalDistribution[μ, σ], x]

Here is the distribution that best fits our data.

NonlinearFit[normdata, model, x, {μ, σ}]

Here we plot the distribution along with the data points.

Plot[%, {x, -3, 3}, Epilog -> Map[Point, normdata]];

Here we look at the representation of the binomial distribution, for which most of the results are well known. BinomialDistribution[n, p] represents the binomial distribution for n trials with probability of “success” p.

Domain[BinomialDistribution[n, p]]

The p.d.f. is expressed in terms of the built-in function Binomial.

binopdf = PDF[BinomialDistribution[n, p], m]

Here we substitute numeric values for m and n.

binopdf /. {m -> 1, n -> 3}

Here are several functions obtained from a binomial distribution with symbolic parameters n and p.

CDF[BinomialDistribution[n, p], m]

Mean[BinomialDistribution[n, p]]

Variance[BinomialDistribution[n, p]]

Skewness[BinomialDistribution[n, p]]

CharacteristicFunction[BinomialDistribution[n, p], t]

The form Random[distribution] generates a random number taken from a specified statistical distribution. Here is a table of 100 random numbers taken from a binomial distribution with parameters 10 and 0.5.

binomialdata = Table[Random[BinomialDistribution[10, 0.5]], {100}];

Needs["Graphics`Graphics`"]; 
    Needs["Statistics`DataManipulation`"]

BarChart[Frequencies[binomialdata], PlotRange -> All];

The function RandomArray creates arrays of random numbers that follow a given statistical distribution. Here is another way to create a list of 100 random numbers that follow a binomial distribution.

RandomArray[BinomialDistribution[10, 0.5], {100}];

BarChart[Frequencies[%], PlotRange -> All];

Here we create a 500 × 2 array of numbers that follow a standard normal distribution, which we interpret as a list of 500 ordered pairs.

RandomArray[NormalDistribution[0, 1], {500, 2}];

Here is a scatter plot of the random array.

ListPlot[%, AspectRatio -> Automatic];

Confidence Intervals and Hypothesis Testing

Normally Distributed Data

The following command generates normally distributed data around zero.

normaldata = RandomArray[NormalDistribution[0, 1], {1000}];

Using BinCounts and BarChart we can see that these data appear to be normally distributed.

dataplot = 
      BarChart[BinCounts[normaldata, {-5, 5, 0.5}], 
        PlotRange -> All];

Here are the sample mean and variance of our random data.

Mean[normaldata]

Variance[normaldata]

Confidence Intervals

The package Statistics`ConfidenceIntervals` allows us to compute confidence intervals for the mean of the population from which data are drawn. The functions ChiSquareCI, NormalCI, ConfidenceLevel, StudentTCI, FRatioCI, VarianceCI, MeanCI, VarianceRatioCI, and MeanDifferenceCI are defined in the package. MeanCI[list, options] returns a list {min, max} representing a confidence interval for the population mean, using the entries in list as a sample drawn from the population.

Needs["Statistics`ConfidenceIntervals`"]

meancilist = MeanCI[normaldata]

By default, MeanCI computes a 95% confidence interval. We change the confidence level by changing the value of the option ConfidenceLevel.

Options[MeanCI]

Here we compute a confidence interval for the mean using a higher confidence level. Notice that the interval becomes larger, as expected.

MeanCI[normaldata, ConfidenceLevel -> 0.975]

We can also do this calculation manually. Here we compute the standard error of the sample mean and the (95^th) quantile of the Student t distribution with 999 (the number of points in the data set, minus one) degrees of freedom.

se = StandardErrorOfSampleMean[normaldata]

q = Quantile[StudentTDistribution[999], (1.95/2)]

Here we compute the mean confidence interval directly.

{Mean[normaldata] - q*se, Mean[normaldata] + q*se}

The result is the same as the result from MeanCI.

% == meancilist

Sometimes we need to calculate confidence intervals based on known parameters, rather than on those estimated from the sample. The first example is another way of doing the calculation done by MeanCI, where we have inserted the mean and its standard error, together with the sample size, while the second example shows what happens with the ordinary normal distribution. We obtain ±1.96, as expected.

StudentTCI[mean, se, dof, ConfidenceLevel-> c] returns a list {min, max} representing a confidence interval at confidence level c for the population mean, based on a sample mean, its standard error se, and dof degrees of freedom. This function is used by MeanCI when the population variance is estimated from the sample.

StudentTCI[-0.017821, 0.102998, 99, ConfidenceLevel -> 0.95]

NormalCI[0, 1, ConfidenceLevel -> 0.95]

We may carry out an analogous series of estimates for population variances.

Variance[normaldata]

VarianceCI[normaldata]

Given two samples, we make comparisons between them using the functions MeanDifferenceCI and VarianceRatioCI.

Hypothesis Testing

We may test hypotheses about the underlying population based on a sample. The output of such a test is a p-value, the probability that the sample estimate is true is given that the value of the hypothesized parameter is correct. The tests defined in the package Statistics`HypothesisTests` are MeanTest, MeanDifferenceTest, VarianceRatioTest, and VarianceTest.

Needs["Statistics`HypothesisTests`"]

Let us test the hypothesis that our data come from a distribution with zero mean. MeanTest returns the p-value, defaulting to carrying out a one-sided test, but can optionally perform a two-sided test.

MeanTest[normaldata, 0, TwoSided -> True]

If is therefore very likely that the data came from a distribution with mean zero. We may extract more detailed information by asking for a full report. We ask for the probability that the data came from a population with mean unity, and suggest a significance level for application of the test. Here we reject the hypothesis that the data come from a normal distribution with mean 1.

MeanTest[normaldata, 1, TwoSided -> True, FullReport -> True, 
      SignificanceLevel -> 0.05]

The various tests that are available parallel the confidence interval functions.

VarianceTest[normaldata, 1, TwoSided -> True, 
      FullReport -> True]

MeanDifferenceTest and VarianceRatioTest are used to compare two sample populations.

5. Related Packages

In addition to the tools presented in this course, Wolfram Research offers more specific add-on packages for data analysis. Two of these add-ons are the Time Series and Experimental Data Analyst packages. Each is available directly from Wolfram Research.

Time Series

This package performs time series analysis of both univariate and multivariate autoregressive moving average (ARMA) types of models. Analysis of seasonal ARMA models and autoregressive integrated moving average (ARIMA) models is also included. The best linear predictions and maximum likelihood used in the maximum likelihood estimations are calculated exactly using the innovations algorithm. In addition, a number of data smoothing functions are introduced.

The functions available in the package Time Series fall into a number of categories.

Needs["TimeSeries`TimeSeries`"];

Time Series Models

The commonly used time series models are represented in generic form and themselves do not evaluate to anything. The models include ARModel, MAModel, ARMAModel, ARIMAModel, and SARIMAModel.

\(ma = 
	
      MAModel[{{{0.1, \(-0.2\)}, {\(-0.2\), 0.1}}, {{0.7, \(-1.2\)}, {0.3, 
              \(-0.2\)}}}, {{1, 0}, {0, 1}}]\)

\(ar\  = 
	
      ARModel[{{{0.2, \ \(-0.3\)}, \ {0.2, \ \(-0.7\)}}, 
          \ {{0.1, \ 0.7}, \ {0.3, \ \(-0.8\)}}}, \ {{1, \ 0}, \ {0, \ 1}}]
        \)

arma  =  ARMAModel[{0.5,  -0.2},  {0.7},  0.8]

arima  =  ARIMAModel[2,  {0.3,  0.9},  {1.2},  1]

sarima  =  
      SARIMAModel[{2,  3},  4,  {0.4,  0.1},  {-0.5},  {0.8},  {}, 
         1]

Extracting Parameters from Time Series Models

These functions allow us to extract information from the time series models. They include ARCoefficients, MACoefficients, SeasonalARCoefficients, SeasonalMACoefficients, SeasonalPeriod, DifferencingOrder, and NoiseParameter.

MACoefficients[arma]

DifferencingOrder[sarima]

Random Sequences and Time Series

We can generate random sequences and time series with the functions RandomSequence and TimeSeries.

RandomSequence[0,  1.2,  5]

TimeSeries[arma,  7]

TimeSeries[ma,  5]

Conversions of Time Series Models

The functions ToARModel, ToMAModel, and ToARMAModel allow us to convert between different time series models.

ToMAModel[ar,  3]

ToARMAModel[arima]

Properties of Time Series Models and Observable Data

We can also query the different properties of the time series models and observable data. This is accomplished with the functions StationaryQ, InvertibleQ, CovarianceFunction, CorrelationFunction, PartialCorrelationFunction, and Spectrum. The functions SmoothedSpectrumS and SmoothedSpectrumL are used with data lists.

StationaryQ[arma]

InvertibleQ[arma]

CovarianceFunction[arma,  2]

CorrelationFunction[ar,  5]

PartialCorrelationFunction[ar,  4]

Estimations of ARMA Models

Given a model, various methods exist to fit the appropriately transformed data to it and estimate the parameters. This is done with the functions YuleWalkerEstimate, LevinsonDurbinEstimate, BurgEstimate, InnovationEstimate, LongAREstimate, and HannanRissanenEstimate.

YuleWalkerEstimate[
      TimeSeries[
        ARModel[{{{0.3, -0.1}, {0.4, 0.2}}, {{-0.5, 0.1}, {-0.2, 
                0.7}}}, {{1, 0}, {0, 1}}], 200], 2]

BurgEstimate[TimeSeries[ARModel[{0.5, -0.2}, 1], 200], 3]

InnovationEstimate[TimeSeries[MAModel[{0.5, -0.2}, 1], 200], 4]

Maximum Likelihood Estimations

MLEstimate gives the maximum likelihood estimate of an ARMA type of model by maximizing the exact likelihood function that is calculated using the innovations algorithm. Since finding the exact maximum likelihood estimate is generally slow, a conditional likelihood is often used. ConditionalMLEstimate gives an estimate of an ARMA model by maximizing the conditional likelihood using the Levenberg-Marquardt algorithm.

\(ConditionalMLEstimate[
	
      TimeSeries[
		
        ARModel[{{{0.3, \ \(-0.1\)}, \ {0.4, \ 0.5}}}, 
          \ {{1, \ 0}, \ {0, \ 1}}], \ 200], 
	1]\)

Predicting Time Series

With the functions BestLinearPredictor and IntegratedPredictor, we can forecast future values of a time series. The best linear predictor is defined as the linear combination of observed data points that has the minimum mean-square distance from the true value. For an ARIMA or SARIMA series the prediction for future values can be obtained using IntegratedPredictor.

model = 
      ARMAModel[{{{0.5, -0.2}, {0.3, -0.3}}}, {{{-0.5, 0.2}, {
              0.8, 0.3}}}, {{1, 0}, {0, 1}}];

BestLinearPredictor[TimeSeries[model,  20],   model,  3]

Data Smoothing

A series of raw data is sometimes transformed to a new series of data before it is analyzed. The purpose of this transformation is to smooth out local fluctuations in the raw data. The functions MovingAverage, MovingMedian, LinearFilter, and ExponentialSmoothing provide this data smoothing capability for us.

Needs["TimeSeries`DataSmoothing`"]

data  = 
       {7.,  4.,  5.,  3.,  8.,  5.,  2.,  3.,  6.,  8.};

MovingAverage[data,  2]

MovingAverage[{{a,  b},  {c,  d},  {e,  f},  {g,  h}},  3]

MovingMedian[data,  3]

LinearFilter[data,  {1/2,  1/2}]

ExponentialSmoothing[data1,  0.8,  0.5]

Experimental Data Analyst (EDA)

Experimental Data Analyst (EDA) is a collection of tools and tutorials designed specifically for the needs of the physical scientist and engineer and students of science and engineering. Included are over thirty programs, tutorials, and discussions of the underlying theory of data analysis, at a variety of levels from beginning undergraduate to practicing researcher. The tutorials make use of a variety of real experimental data drawn from a variety of different fields, and the data are also included with the pack.

EDA is divided into eight areas which are discussed at length in the provided documentation. Here is a brief summary and some examples from each section.

Needs["EDA`Master`"];

Data Input

There are a number of methods of getting data into Mathematica so that it may be analyzed. This section discusses techniques to read and write files containing data sets and introduces the EDA functions ImportData and ExportData.

Error Analysis

In this section, EDA supplies functions and constructs to simplify error analysis. The EDA programs discussed are AdjustSignificantFigures, CombineWithError, DivideWithError, PlusWithError, PowerWithError, Quadrature, SubtractWithError, and TimesWithError. Also, EDA defines Data and Datum constructs for doing propagation of errors.

Usually, errors of precision are probabilistic. This means that the experimenter is saying that the actual value of some parameter is probably within the specified range. For example, if the error is the standard deviation then the probability is about 68%; if the error is twice the standard deviation the probability is 95%, etc. If we have two variables, say x and y, and wants to combine them, we want the error in the combination to preserve this probability. The correct procedure to do this is to combine errors in quadrature (i.e. the square root of the sum of the squares).

Quadrature[a,  b,  c]

EDA includes programs to combine data including errors using quadrature. These are the *WithError functions.

p = {{70, 0.04}, {72, 0.04}, {74, 0.04}, 
	    {76, 0.04}, {78, 0.04}, {80, 0.04}, 
	    {82, 0.04}, {84, 0.04}, {86, 0.04}, 
	    {88, 0.04}, {90, 0.04}};

v = {{11.28156820762763, 0.031}, 
     {10.98366389613071, 0.031}, 
     {10.6916855829648, 0.031}, 
     {10.4292141352925, 0.031}, 
     {10.15135483487289, 0.031}, 
     {9.9030178446771, 0.031}, 
     {9.69180367127115, 0.031}, 
     {9.45284578420525, 0.031}, 
     {9.23770984228389, 0.031}, 
     {9.03759086299312, 0.031}, 
     {8.84919072382683, 0.031}};

TimesWithError[p, v]

PowerWithError[p, 2]

Fitting Data to Linear Models (Straight Lines and Curves)

This section discusses what a linear fit is, least square techniques, evaluation of the quality of the fit, etc.; some familiarity with this material is assumed in all later sections. It also discusses the EDA program LinearFit that performs linear least square fits to data; included is discussion of the related EDA programs ShowLinearFit and ToLinearFunction.

LoadData[Ganglion]

result = LinearFit[GanglionData, {0, 1}, a]

Fitting Data to Arbitrary Models

This section introduces the EDA program FindFit and the related programs ShowFitFunction and ToFitFunction. Also discussed are some convenience programs that define peak shapes; these are: BreitWigner, Galatry, Gaussian, Lorentzian, PearsonVII, RelativisticBreitWigner, and Voigt.

LoadData[Cobalt60]

\(FindFit[Cobalt60Data, 
	
      a[0] + a[1]\ chan + Gaussian[chan, ampl, x0, sigma], 
      chan, {{a[0], 234}, {a[1], \(-0.11\)}, {ampl, 140}, {x0, 1830}, {sigma, 
          45}}]\)

Needs["Graphics`Legend`"]; Plot[{Gaussian[x, 10, 10, 2.5], BreitWigner[x, 80, 10, 5]}, {x, 0, 
        20}, PlotStyle -> {GrayLevel[0], Dashing[{0.03}]}, 
      PlotLegend -> {"Gaussian", "BreitWigner"}, 
      LegendSize -> {1, 0.4}];

Eliminating Noise in Data and Filling in Missing Values

This section discusses techniques to eliminate noise in data and to fill in missing values. It also discusses the related topic of fitting data when a model is not available or appropriate. The chapter introduces the EDA programs SmoothData, LoessFit, and FillData. There is also a tutorial on using Mathematica built-in programs to smooth data.

Fitting Data with “Wild” Data Points

This section discusses techniques to fit data to lines and curves when one or more of the data points may be “wild”; least square techniques should not be used in this case since a wild data point can seize control of the least square fit. The chapter introduces the EDA programs RobustCurveFit and RobustLineFit.

RobustLineFit[GanglionData, a]

SeedRandom[1234]; mydata = Table[1.5 + Random[Real, {-0.1, 0.1}], {20}]; mydata = Insert[mydata, 100, 5];

LinearFit[mydata, {0}, a, ReturnErrors -> False]

RobustCurveFit[mydata, {0}, a]

What to Do with Unknown Relations in a Data Set

This section discusses what to do when the relations in a data set are not known. Graphical techniques are explored and  the EDA program EDAListPlot is introduced. The EDA programs Histogram and BoxPlot are also discussed.

mydata1 = Table[{x, Sin[x]}, {x, 0, π, (π/20)}]; mydata2 = Table[{x, {Cos[x], 0.2}}, {x, 0, π, (π/25)}];

EDAListPlot[mydata1, mydata2];

LoadData[KPiPi]

Histogram[(Sqrt[Last])/@KPiPiData, 75, 0.5, 1.25, 
      PlotLabel -> "(Pi+, Pi-) mass", ShowAll -> True];

LoadData[Darwin]

BoxPlot[{CrossFertilizedData, Take[CrossFertilizedData, -12], 
        SelfFertilizedData, Take[SelfFertilizedData, -12]}, 
      Prolog -> PointSize[0.02], 
      FrameTicks 
        -> {{{1, "Cross"}, {2, "Part Cross"}, {3, "Self"}, {
              4, "Part Self"}}, Automatic}];