Chapter 2

Descriptive Statistics

The X-Files data

Mathematica has functions to perform many statistical operations, but they are not generally loaded when Mathematica starts.  Rather related functions are bundled into packages which can be loaded as a group before they are needed.  Measures of center and spread are in a package called `DescriptiveStatistics` and it can be loaded with the following command.

Needs["Statistics`DescriptiveStatistics`"];

Once loaded the descriptive statistic fuctions are available. Consider the data

season1
       = {7.4, 6., 6.8, 5.9, 6.2, 5.6, 5.6, 6.2, 6.1, 5.1, 6.4, 6.4, 6.2, 
        6.8, 7.2, 6.8, 5.8, 7.1, 7.2, 7.5, 8.1, 7.7, 7.4, 8.3}

Mean[season1]

Median[season1]

Mode[season1]

Variance[season1]

StandardDeviation[season1]

SampleRange[season1]

InterquartileRange[season1]

Quantile[season1, .5]

Quantile[season1,  .25]

Quantile[season1,  .75]

To make bar charts of the data, the packages Statistics`DataManipulation` and Graphics`Graphics` are needed.

Needs["Statistics`DataManipulation`"]

The Frequencies command counts the number of occurences of each piece of data.

freq = Frequencies[season1]

We complete the histogram by making a bar chart of the frequencies.  This is a slightly different method of histogram construction than mentioned in the text.  There the bars covered the intervals of interest; here the bars are centered on the midpoints of the intervals and the bar widths have no meaning.  Other packages may use still other methods.

Needs["Graphics`Graphics`"]

BarChart[freq];

The BinCounts function can be used to bin the data.  Here we group the data into 7 bins.

freq = BinCounts[season1, {5.09, 8.30, (8.31 - 5.09)/7}]

We also need the midpoint of each interval.

midpoints = 
      Table[5.0 - ((8.3 - 5.1)/(2*7)) + 
          n (8.3 - 5.1)/(7), {n, 1, 7}]

BarChart[Transpose[{freq, midpoints}]];

Next we look at constucting box plots of data.  Mathematica does not come with a boxplot function so we create on using graphics elements.  Note that we use Statistics`DescriptiveStatistics` package.  If it has been already loaded, it need not be reloaded.

season1
       = {7.4, 6., 6.8, 5.9, 6.2, 5.6, 5.6, 6.2, 6.1, 5.1, 6.4, 6.4, 6.2, 
        6.8, 7.2, 6.8, 5.8, 7.1, 7.2, 7.5, 8.1, 7.7, 7.4, 8.3}

Needs["Statistics`DescriptiveStatistics`"];

\(SingleBoxplot[data_, center_, \ width_] := 
      \((
		Clear[a, b, c, d, e, f, g]; 
		a = Min[data]; 
        b = \(Quartiles[data]\)[[1]]; 
        c = \(Quartiles[data]\)[[2]]; 
        d = \(Quartiles[data]\)[[3]]; 
        e = Max[data]; f = center - width; g = center + width; 
		
        Graphics[{
					Line[{{f, e}, {g, e}}], 
					
            Line[{{f, d}, {g, d}}], 
					Line[{{f, c}, {g, c}}], 

            					Line[{{f, b}, {g, b}}], 
					
            Line[{{f, a}, {g, a}}], 
					Line[{{f, b}, {f, d}}], 

            					Line[{{g, b}, {g, d}}], 
					
            Line[{{center, a}, {center, e}}]}])\)\)

\(\(s1 = SingleBoxplot[season1, 1, .25]; 


    Show[s1, Axes -> True, Ticks -> {{1}, Automatic}, 
      AxesOrigin -> {.5, Automatic}]; 
	\)\)

season2
       = {9.8, 9.3, 8.7, 8.2, 8.5, 9.2, 9, 9.1, 8.6, 9.9, 8.5, 9.7, 8.8, 
        10.2, 10.8, 9.8, 10.7, 9.6, 10.2, 9.8, 7.9, 8.5, 8.1, 9, 9.6}

\(\(s2 = SingleBoxplot[season2, 1, .25]; 


    Show[s2, Axes -> True, Ticks -> {{1}, Automatic}, 
      AxesOrigin -> {.5, Automatic}]; 
	\)\)

SingleBoxPlot can be used to make side-by-side box plots.

\(\(s1 = SingleBoxplot[season1, 1, .25]; 

    s2 = SingleBoxplot[season2, 2, .25]; 

    Show[s1, s2, Axes -> True, Ticks -> {{1, 2}, Automatic}, 
      AxesOrigin -> {.5, Automatic}]; 
	\)\)

Stochastic Simulations

The Sandhill Cranes

Random numbers can be generated from any standard distribution and used in simulations.  Before these are used the package Statistics`ContinuousDistributions` must be loaded.

Needs["Statistics`ContinuousDistributions`"];

We simulate the crane population with demographic stochasticity.  Here the birth rate is normal with mean 0.5 and standard deviation 0.03 and the death rate is normal with mean 0.1 and standard deviation 0.08.

Here is one simulation.

\(Clear[n, t]; 
n[0] = 100; 

    n[t_] := \(n[t] = 
        \((Random[NormalDistribution[0.5, 0.03]] - 
				
              Random[NormalDistribution[0.1, 0.08]] + 1)\)\ n[t - 1]\); 
    ListPlot[Table[{t, n[t]}, {t, 0, 10}], PlotJoined -> True]; \)

Next we create a table of simulations and show them on one plot.

Clear[s]; \(s[i_] := 
      \((Clear[n, t]; n[0] = 100; 
        n[t_] := 
          \(n[t] = 
            \((Random[NormalDistribution[0.5, 0.03]] - 
				
                  Random[NormalDistribution[0.1, 0.08]] + 1)\) n[t - 1]\); 
        ListPlot[Table[{t, n[t]}, {t, 0, 10}], PlotJoined -> True, 
          DisplayFunction -> Identity])\); 

    Show[Table[s[i], {i, 1, 100}], DisplayFunction -> $DisplayFunction, 
      AxesLabel -> {"Years", "Population"}]; \)

We will not replicate the environmental stochasticity model, but indicate that uniform random numbers from [a,b] can be generated with Random[Real, {a, b}]

Thus to generate random numbers from the unit interval we use.

Random[Real, {0, 1}]

The model can be completed with an If statement.  The syntax is If[test, trueresult, falseresult].
The following will return .3 forty percent of the time and .5 the remaining sixty percent.

If[Random[Real, {0, 1}] < .4,  .3,  .5]