Chapter 2

Descriptive Statistics

Needs["Statistics`DescriptiveStatistics`"];

house
       = {70, 65, 80, 91, 105, 143, 150, 145, 120, 92, 100, 95, 83, 139, 510, 
        101, 675, 84, 580}ast = {19, 25, 22, 23, 24, 19, 22, 19, 24, 24, 27, 18, 20, 23, 25, 28, 
        10, 18, 23, 19, 21, 17, 23, 26, 24, 26, 21, 14, 27, 18, 19, 24, 23, 
        21, 23, 20, 22, 25, 25, 16, 24, 19, 22, 15, 22, 26, 19, 8, 22, 24}gpigs
       = {40,  45,  53,  56,  56,  57, 58, 62, 64, 73, 74, 79, 78, 80, 
        81, 81, 81, 82, 83, 83, 84, 85, 89, 89, 91, 92, 92, 97, 99, 99, 99, 
        99, 101, 102, 102, 102, 103, 104, 107, 108, 109, 113, 114, 118, 121, 
        123, 126, 128, 137, 138, 134, 144, 145, 147, 156, 162, 174, 174, 179, 
        184, 191, 198, 208, 214, 247, 249, 328, 383, 403, 511}cereal
       = {1622, 1579, 1603, 1576, 1589, 1625, 1588, 1589, 1562, 1595, 1589, 
          1599, 1603, 1594, 1629, 1606, 1617, 1597, 1616, 1614, 1600, 1616, 
          1596, 1589, 1608, 1594, 1596, 1608, 1606, 1617, 1618, 1610, 1598, 
          1617, 1593, 1587, 1592, 1603, 1605, 1598, 1599, 1592, 1622, 1625, 
          1564, 1603, 1605, 1626, 1614}/100.00

The house data: #1, p. 90

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.

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

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

Mean[house]*1.0Median[house]\(Clear[n]; 
Nbin = 15; 
smallBit = .01\)minMin = Min[house] - smallBitmaxMax = Max[house] + smallBitfreq = 
      BinCounts[house, {minMin, Max[house], (maxMax - minMin)/Nbin}]midpoints = 
      Table[minMin + 
          (2  n - 1) (Max[house] - Min[house])/(2*Nbin), {n, 1, 
          Nbin}]BarChart[Transpose[{freq, midpoints}], Axes -> Automatic, 
      PlotRange -> {0,  12}]

Other data sets: #3, p. 90

\(Clear[n]; 
Nbin = 15; 
smallBit = .01\)minMin = Min[gpigs] - smallBitmaxMax = Max[gpigs] + smallBitfreq = 
      BinCounts[gpigs, {minMin, Max[gpigs], (maxMax - minMin)/Nbin}]midpoints = 
      Table[minMin + 
          (2  n - 1) (Max[gpigs] - Min[gpigs])/(2*Nbin), {n, 1, 
          Nbin}]BarChart[Transpose[{freq, midpoints}], Axes -> Automatic, 
      PlotRange -> {0,  30}]

\(Clear[n]; 
Nbin = 7; 
smallBit = .001\)minMin = Min[cereal] - smallBitmaxMax = Max[cereal] + smallBitfreq = 
      BinCounts[cereal, {minMin, Max[cereal], (maxMax - minMin)/Nbin}]midpoints = 
      Table[minMin + 
          (2  n - 1) (Max[cereal] - Min[cereal])/(2*Nbin), {n, 1, 
          Nbin}]BarChart[Transpose[{freq, midpoints}]]

\(Clear[n, list]; 
Nbin = 7; 

    season1 = {7.4, 6.9, 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}\)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}\(list = season2; 
smallBit = .001\)minMin = Min[list] - smallBitmaxMax = Max[list] + smallBitfreq = 
      BinCounts[list, {minMin, Max[list], (maxMax - minMin)/Nbin}]midpoints = 
      Table[minMin + 
          (2  n - 1) (Max[list] - Min[list])/(2*Nbin), {n, 1, 
          Nbin}]BarChart[Transpose[{freq, midpoints}]]

\(Clear[n, list]; 
Nbin = 7; 
list = Join[season2, season1]; 

    smallBit = .001\)minMin = Min[list] - smallBitmaxMax = Max[list] + smallBitfreq = 
      BinCounts[list, {minMin, Max[list], (maxMax - minMin)/Nbin}]midpoints = 
      Table[minMin + 
          (2  n - 1) (Max[list] - Min[list])/(2*Nbin), {n, 1, 
          Nbin}]BarChart[Transpose[{freq, midpoints}]]

Xfiles data: #4, p. 91

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 Nbin bins.

\(Nbin = 7; 
smallBit = .01\)minMin = Min[season1] - smallBitmaxMax = Max[season1] + smallBitfreq = 
      BinCounts[season1, {minMin, Max[season1], (maxMax - minMin)/Nbin}]

We also need the midpoint of each interval.

midpoints = 
      Table[minMin - (Max[season1] - Min[season1])/(2*Nbin) + 
          n (Max[season1] - Min[season1])/Nbin, {n, 1, Nbin}]

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.

\(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}]; 
	\)\)

\(s2 = SingleBoxplot[season2, 2, .25]; 
Min[season2]\)Quantile[season2, .25]Quantile[season2, .5]Median[season2]Quantile[season2, .75]Max[season2]

\(
\)

\(\(\(Show[s2, Axes -> True, Ticks -> {{2}, Automatic}, 
      AxesOrigin -> {1.5, Automatic}]; \)
	\)\)

\(
\)

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

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

Mean and Variance of a random variable: #6, p. 91

Can't be done simply in Mathematica, so I'll just write this one out. The idea is that
E(aX+bY)=∫∫ (ax+by)p(x,y)dx dy
where p(x,y) is the joint density for the random variables X and Y. This is equal to
∫∫ axp(x,y)dx dy + ∫∫ byp(x,y)dx dy
and if we do these integrals in the appropriate order, then things become easy:
E(aX+bY)=∫ ax (∫ p(x,y)dy) dx  + ∫ by (∫ p(x,y)dx ) dy
		=∫ ax f(x) dx  + ∫ by g(y) dy
where f(x) is the pdf of X, and g(y) is the pdf of Y. Thus, 
E(aX+bY)	= E(aX)  + E(aY)

The uniform distribution: #7, p. 91

Let's do a few problems using the uniform: first we need to load in the continuous distributions:

<< Statistics`ContinuousDistributions`

udist  =  UniformDistribution[alpha, beta]pdf  =  PDF[udist,  x]

As for the mean, this is the weighted integral of the pdf times x over all the reals:

Simplify[Mean[udist]]Simplify[Integrate[x/(beta - alpha), {x, alpha, beta}]]

Now for the variance:

Simplify[Variance[udist]]Simplify[
      Integrate[
        (x - ((alpha + beta)/2))^2/(beta - alpha), {x, alpha, beta}]]

The normal distribution: #8, p. 91

ndist  =  NormalDistribution[mu, sigma]pdf  =  PDF[ndist,  x]

The mean is the weighted integral of the pdf times x over all the reals. You can do this integral by parts.

Simplify[Mean[ndist]]Simplify[Integrate[x* pdf, {x, -∞, ∞}]]

The variance is the weighted integral of the pdf times (x-mu)^2 over all the reals. You can do this integral by parts.

Simplify[Variance[ndist]]Simplify[
      Integrate[(x - mu)^2*pdf, {x, -∞, ∞}]]

To find max/mins, we differentiate the pdf and solve for pdf=0:

Simplify[Dt[pdf, x]]

Since the mu and sigma are constants, their derivatives are zero, and the only zero of this function is x=mu.

Simplify[Dt[pdf, {x, 2}]]

Again, since the mu and sigma are constants, most of these terms are zero. Here's how you indicate that mu and sigma are constants to Mathematica:

Dt[pdf,x, Constants->{mu, sigma}]Simplify[Dt[pdf,{x, 2}, Constants->{mu, sigma}]]

Now let's use Solve to find those zeros:

Solve[Dt[pdf,x, Constants->{mu, sigma}]==0,x]Solve[Simplify[Dt[pdf,{x, 2}, Constants->{mu, sigma}]]==0,x]

As for the 68-95-99 rule, we simply need to integrate the pdf between the required intervals:

Integrate[pdf, {x, mu - sigma, mu + sigma}]*1.0Integrate[pdf, {x, mu - 2  sigma, mu + 2  sigma}]*1.0Integrate[pdf, {x, mu - 3  sigma, mu + 3  sigma}]*1.0

Guinea Pig problem

For this problem I asked you to test the hypothesis that the guinea pig data is normal using a chi-square test. For that, we create some bins assuming a normal distribution, and compare observed to expected:

\(dist\  = \ NormalDistribution[Mean[gpigs], StandardDeviation[gpigs]]; 

    
nBins = 5; 
step = 1.0/nBins; \)

cutpoints is the location of the first 100/nBins% of the data, 
the second 100/nBins% of the data, etc., from the theoretical distribution:

cutpoints = Quantile[dist, Table[n*step,  {n, 1, nBins - 1}]]

Now RangeCounts tells us how many observations fall in each bin:

chifreqs = RangeCounts[gpigs, cutpoints]

Because of the way we've chosen our bins, the same number are expected in each bin:

expected = Length[gpigs]/nBins*1.0

Note that the chi-square statistic is actually just the sum of the chisq observeds minus  expecteds:

chisq = Mean[(chifreqs - expected)^2/expected]*nBins

And the p-value of the statistic is given by the area under the curve from chisq to infinity (which is the difference below):

1 - CDF[ChiSquareDistribution[nBins - 1], chisq]

Hence we can safely reject the hypothesis of a normal distribution for this variable. It's clearly unlikely!