Mathematica statistics tips:

Mathematica needs to load in the requisite statistical functions, which we do as  follows:

<< Statistics`ContinuousDistributions`<< Statistics`DataManipulation`<< Graphics`Graphics`

Once loaded, we can access the distributions of all kinds. We're interested 
in the normal, particularly:

PDF[NormalDistribution[mu,  sigma],  x]\(dist\  = \ NormalDistribution[0, 1]; \ \ \ \ \ 
     (*\ this\ is\ the\ normal\ with\ mu = 0, \ sigma = 1\ *) 

    pdf\  = \ PDF[dist, \ x]\)cdf  =  CDF[dist,  x]

Have a look at the plots, which should be familiar by now:

Plot[pdf, {x,  -3,  3}]Plot[cdf, {x,  -3,  3}]Plot[{pdf, cdf}, {x,  -3,  3}]

But we can play with lots of other distributions, if we need to:

PDF[GammaDistribution[alpha,  lambda],  x]Plot[PDF[GammaDistribution[1.9,  1.9],  x], {x,  0, 15}]

The quantile function is the inverse of the cdf function:

Quantile[dist,  CDF[dist,  .68]*1.0]    
      \(CDF[dist, Quantile[dist, \ .68]]\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
     (*\ should\ return\ .68\ *) 
\)\(quant\  = \ Quantile[dist, x]; 

    Plot[{cdf, quant, x}, {x, \ \(-3\), \ 3}]\)Plot[quant, {x,  -1,  1}]

We can compute Means and Standard Deviations of distributions, as well as 
data:

Mean[dist]StandardDeviation[dist]

As we discussed in class, it's often necessary to generate random  variables. We've seen that it's possible to generate random variables for any  given distribution (pdf), if we have a uniform random variable generator and  the cdf. Let's compare the normal random number generator of Mathematica and "our" random number generator:

\(Random[dist]\ \ \ \ \ \ \ \ \ 
     (*\ the\ function\ Random\ produces\ random\ numbers\ from\ a\ given\ 
        distribution\ *) 
\)\(numRands = 1000; 

    theirRandomNormals\  = \ Table[Random[dist], \ {numRands}]; 

     (*\ this\ is\ a\ list\ of\ 1000\ random\ numbers\ *) 


     (*\ new\ we' ll\ generate\ our\ random\ normals\ using\ a\ uniform\ 
        distribution, \ 
      and\ the\ inverse\ of\ the\ cdf\ of\ the\ standard\ \(normal : \)\ *) 

    unifo\  = \ UniformDistribution[0, \ 1]; 

    ourRandomNormals\  = \ 
      Quantile[dist, \ \ Table[Random[unifo], \ {numRands}]]*1.0; \)

What should we expect to see here?

Mean[theirRandomNormals]Mean[ourRandomNormals]StandardDeviation[theirRandomNormals]StandardDeviation[ourRandomNormals]

Let's look at a couple of histograms of the data, so that we can compare  the two. Mathematica doesn't have a nice histogram plot type, but we can kludge one together  this way:

\(\  (*\ we' re\ looking\ at\ binNum\ bins, \ 
      from\ binMin\ to\ binMax\ *) 
binMin = \(-4\); 
binMax = 4; 

    binNum = 8; 
binSize = \((binMax - binMin)\)/binNum; 


     (*\ we\ now\ compute\ the\ frequencies\ in\ the\ various\ 
        \(bins : \)\ *) 

    ourFreq\  = \ 
      BinCounts[ourRandomNormals, \ {binMin, \ binMax, \ binSize}]; 

    theirFreq\  = \ 
      BinCounts[theirRandomNormals, \ {binMin, \ binMax, \ binSize}]; 


    midpoints\  = \ Table[binMin + n*binSize, \ {n, 0, binNum - 1}]; 


    BarChart[Transpose[{ourFreq, \ midpoints}]]\)BarChart[Transpose[{theirFreq,  midpoints}]]

Let's look at two plots:
	a plot of the sample joint density of the two, and
	a plot of the quantiles of the two samples:

ListPlot[Transpose[{ourRandomNormals, theirRandomNormals}]]ListPlot[Transpose[{Sort[ourRandomNormals], Sort[theirRandomNormals]}]]

We compute cutpoints for dividing up the data into our equiprobable bins:

\( (*\ 
      cutpoints\ is\ the\ location\ of\ the\ first\ 10  %\ of\ the\ data, \ 
      20  %\ of\ the\ data, \ \(etc.\), \ 
      from\ the\ theoretical\ \(distribution : \)\ *) 
step = .10; 

    nBins = 10; 

     (*\ cutpoints = 
        Quantile[dist, {.1, .2, \ .3, \ .4, .5, .6, .7, .8, .9}]\ *) 

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


     (*\ here' s\ how\ you\ can\ access\ a\ particular\ element\ *) \)\(cutpoints[\([1]\)]


     (*\ now\ RangeCounts\ tells\ us\ how\ many\ fall\ in\ each\ 
        \(bin : \)\ *) \)chifreqs = RangeCounts[theirRandomNormals, cutpoints]

We can now compute chi-square approximations for the data:

\(expected = numRands/nBins*1.0; 


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

     (*\ note\ that\ this\ is\ actually\ just\ the\ sum\ of\ the\ chisq\ 
          \((observeds\ minus\ expecteds)\)^2/expecteds\  - \ 
        I\ just\ couldn' t\ figure\ out\ how\ to\ use\ the\ Sum\ 
          \(function!\)\ *) 
\)chipdf = PDF[ChiSquareDistribution[9],  x]chicdf = CDF[ChiSquareDistribution[9],  x]1 - CDF[ChiSquareDistribution[9], chisq]Plot[chipdf, {x, 0 ,  3*chisq}]Plot[{chipdf, chicdf}, {x, 0 ,  3*chisq}]

Now here's the data for Project 2.3, p. 95:

\(recruitment
       = {\(-47.8\), 70.3, 100.6, 35.3, 62.0, 89.3, 34.3, 57.8, 13.0, 112.8, 
        239.4, \(-9.3\), 61.7, 99.6, 
57.9, 49.2, 48.9, 61.1, 73.4, 43.6, 
        87.7, 49.2, 72.1}
\)\(survival
       = {0.90, 0.64, 0.85, 0.57, 0.58, 0.64, 0.59, 0.51, 0.69, 0.49, 0.61, 
        0.54, 0.50, 0.43, 0.75, 
0.54, 0.50, 0.50, 0.70, 0.54, 0.48, 0.63, 
        0.55, 0.47}
\)population  = 
       {273.8, 126.8, 178.1, 201.9, 151.4, 159.0, 183.5, 126.5, 145.2, 85.0, 
        164.3, 328.1, 155.7, 127.4,  194.2, 161.2, 129.0, 112.8, 137.9, 
        147.0, 113.2, 157.5, 135.7, 134.2}