Using Mathematica for statistics.

D. E. Soper and S. James Remington
28 March 2005

Built in functions

Mathematica has functions that calculate something,

Sin[0.1]

Mathematica has functions that do something,

Plot[Sin[x],{x,0,3}]

The functions that do something generally also return a "value"

plot1 = Plot[Sin[x],{x,0,2 Pi}]

plot2 = Plot[Cos[x],{x,0,2 Pi}]

Show[plot1,plot2]

To find out about Mathematica functions (or anything else in Mathematica), try the Help browser.

Lists

Mathematica likes to work with lists, where a list is a bunch of things enclosed  in braces and separated by commas. For example a list of numbers could be {1,7,4,4,48}. Here, we make a list of the form  {{0,0},{1,0},{2,0),(3,0},...} with the function "Table." Each list element is itself a list of the  form {n,0}.

results = Table[{n-1,0},{n,1,10}]

In listspeak, the first element of the 5th element of this list is

results[[5,1]]

while the second element of the 5th element of this list is

results[[5,2]]

To make a [[  you can type "escape", "[","[", "escape",while to make a ]] you can type "escape", 
"]","]", "escape". You can also use

results[[5,1]]

but it's not as pretty.

Write your own functions

You can write your own function. The variable name, here x, is denoted x_ 
indicating that x can take any value. You collect variables internal to 
the function, in this case y, as the first part of the Module.

xplusone[x_] := Module[{y},y = x + 1]

Try it.

xplusone[4]

xplusone[t]

A function can have several steps, separated by semicolons. The "answer" 
is the result of the last step.

ff[r_] := 
Module[{x,y,z},
x = r + 1;
(* This is a comment *)
y = 1/x + 5;
(* You can put a comment anywhere *)
z = Sin[y]
]

q = ff[t]

ff[0.1]

A function can do something. The "answer" is the result of the last 
step.

fff[r_] := 
Module[{x,y,z},
Print["Starting with r = ",r];
x = r + 1;
Print["Add 1: ",x];
y = 1/x + 5;
Print["Invert and add 5: ",y];
z = Sin[y];
Print["Take the sine of that: ",z];
z
]

fff[t]

Loops

The instruction Do[ something, {n,1,8}] does what ever "something" is 
with n = 1, then does it again with n = 2, then does it again with n = 3, et cetera. It stops after it does "something" with n = 8. (I gave this 
Module an empty list of internal variables.)

sorry[sin_,nmax_] :=
Module[{},
  Do[
    Print[n, ". I will not ",sin," any more."],
    {n,1,nmax}
    ]
]

sorry["write on the wall",3]

We could use a loop to calculate n!.

myfactorial[n_] :=
Module[{x},
  x = 1;
  Do[
  x = x * j,
  {j,1,n}];
  x
  ]

myfactorial[4]

Actually, Mathematica already knows about this function..

Factorial[4]

4!

Logical functions

Some functions give values True or False.

4 < 5

4 > 5

Branches

You can use logical functions and the If[ ...] construct to get the 
computer to do different things depending on the outcome of some test. This instruction will set z to 3 if n is less than 4 and will set z to 0 otherwise.

n = 2

If[n<4,z=3,(*else*)z=0]

If you leave out the "else" part, it does nothing when the condition is 
not satisfied.

n = 3

If[n<4,Print["n<4"]]

We can use this in a function designed to do one thing if x < 0.3 and 
another otherwise.

h[x_] :=
Module[{p},
  p = 0.3;
  If[x<p,
    Print[x," is less than ",p,"."],
    (*else*)
    Print[x," is greater than or equal to ",p,"."]
    ]
    ]

h[0.2]

Random numbers

You can generate a random number between 0 and 1. Each time you do this, 
the computer chooses a different random number. However, the sequence of 
numbers is always the same when you start from "scratch".

x = Random[]

x  =  Random[]

You can do something at random

You can use the function random to take one action with probability p  that you specify and something else with probability 1-p.

takeachance[p_] :=
Module[{r},
r = Random[];
If[r<p,
   Print["Something good happened.  ",r],
   (* else *)
   Print["Something bad happened.   ", r]
  ]
]

takeachance[0.7]

You can do a little statistical simulation

You can use the function random to take one action with probability p  that you specify and something else with probability 1-p.

chance[p_] :=
Module[{r,x},
r = Random[];
If[r<p,
   x=True, (* x = True with probability p *)
   (* else *)
   x=False (* x = False with probability 1-p *)
  ];
x (* Return the value of x *)
]

chance[0.7]

numberT[p_,n_] :=
Module[{sumofT},
sumofT = 0;
Do[(*j from 1 to n*)
  x = chance[p];
  If[x,sumofT = sumofT + 1],
  {j,1,n}
  ];
sumofT
]

numberT[0.3,10]

This is a lot easier than flipping coins! We find that numberT[p,n] is typically something like p*n, but there is some variation. The variation 
as a fraction of p*n is smaller when n is large.

You can make a graph

Let us set some parameters. (You can always change them.)

p = 0.3;

ntot = 10;

ntries = 300;

Now we do a simulation by flipping a (biased) coin ntries times. Each flip gives T or F and we record how the fraction of the  time that we get exactly n results T. We create a an empty list "results" to hold the results,  in the form

     { n, number of trials in which we got T n times/ total number of trials}.

results = Table[{n-1,0},{n,1,ntot+1}]

Do[
y = numberT[p,ntot];
results[[y+1,2]] = results[[y+1,2]] + 1/ntries,
{j,1,ntries}]

results

Check that our fractions add to one.

Sum[results[[n,2]],{n,1,Length[results]}]

We can plot our results. ( I plot from 0 to 10, but you can change that.)

plot1 = ListPlot[results,
PlotRange -> {{0,10},All},
PlotStyle -> PointSize[0.02]]

Here is the theoretical probability for getting T n times. (You should be 
able to derive this.)

prob[n_] := ntot!/n!/(ntot-n)! * p^n * (1-p)^(ntot-n)

Here is a graph of it.

plot2 = Plot[prob[n], {n,0,10}]

We can overlay the plots. In the simulation, the fraction of times that we 
get the result n is not the probability predicted because we didn't do the 
experiment many times. But if we chose a big ntries, the two plots should 
be close.

Show[plot1,plot2]