Using Mathematica 
                                    in 
                 Probability and Statistics

John Kinney

Rose-Hulman Institute of Technology

Statistics - Single Sample

The data gives the weights (in ounces) of zinc coatings of 75 galvanized iron sheets (Guttman, Wilks, and Hunter, Introductory Engineering Statistics, John Wiley  & Sons,1982)

zincdata
       = {1.47, 1.62, 1.52, 1.77, 1.55, 1.53, 1.38, 1.37, 1.64, 1.46, 1.63, 
        1.54, 1.53, 1.34,  1.60,  1.64,  1.60,  1.60, 1.38,  1.73,  
        1.70,  1.60,  1.60,  1.48,  1.51,  1.53,  1.59,  1.50,  1.60, 
         1.54,  1.57,  1.51,  1.58,  1.58,  1.32,  1.62,  1.47,  
        1.42,  1.45,  1.34,  1.44,  1.56,  1.48,  1.55,  1.64,  1.57, 
         1.44,  1.56,  1.39,  1.65,  1.62,  1.53,  1.47,  1.34,  
        1.58,  1.49,  1.56,  1.54,  1.57,  1.67,  1.47,  1.63,  1.49, 
         1.44,  1.35,  1.53,  1.38,  1.46,  1.44,  1.47,  1.43,  
        1.64,  1.50,  1.61,  1.42,  1.57,  1.75,  1.47,  1.64};

<< Statistics`

LocationReport[zincdata]

DispersionReport[zincdata]

Quartiles[zincdata]

Other statistics available include GeometricMean,  TrimmedMean, Quantiles, VarianceOfSampleMean, QuartileDeviation, and many others.

<< Graphics`Graphics`

Min[zincdata]

Max[zincdata]

freq = BinCounts[zincdata, {Min[zincdata], Max[zincdata], 0.05}]

tix = Table[1.32 + k*(0.05), {k, 0, 8}]

BarChart[Transpose[{freq, tix}], 
      AxesLabel -> {"Thickness", "Frequency"}]

Confidence Intervals

MeanCI[zincdata]

MeanCI[zincdata, ConfidenceLevel -> 0.935, KnownVariance -> 0.02]

VarianceCI[zincdata]

VarianceCI[zincdata, ConfidenceLevel -> 0.9475]

Hypothesis Testing

MeanTest[zincdata, 1.56]

MeanTest[zincdata, 1.56, KnownVariance -> 0.015]

VarianceTest[zincdata, 0.015]

Statistics - Two Samples

zincdata2
       = {1.49, 1.77, 1.34, 1.89, 1.56, 1.74, 1.78, 1.67, 1.84, 1.22, 1.76, 
        1.73, 1.85, 1.94, 1.56}

LocationReport[zincdata2]

DispersionReport[zincdata2]

MeanDifferenceCI[zincdata, zincdata2]

MeanDifferenceCI[zincdata, zincdata2, EqualVariances -> True]

MeanDifferenceCI[zincdata, zincdata2, Known Variance -> {0.02, 0.04}, 
      ConfidenceLevel -> 0.99]

VarianceRatioCI[zincdata, zincdata2]

MeanDifferenceTest[zincdata, zincdata2, 0]

Statistics - Regression

The data represent the number of surviving bacteria 
      (in hundreds) as estimated by plate counts in an experiment with marine 
        bacteriumfollowing exposure to 200 - kilovolt X - 
        rays for periods ranging from t =1 to 15 in intervals of 6 minutes.The data comes from Regression   Analysis By Example by Chatterjee and Price,  
    John Wiley & Sons,  1991.

bactdata
       = {{1, 355}, {2, 211}, {3, 197}, {4, 166}, {5, 142}, {6, 106}, {7, 
          104}, {8, 60}, {9, 56}, {10, 38}, {11, 36}, {12, 32}, {13, 21}, {
          14, 19}, {15, 15}}

ListPlot[bactdata, AxesLabel -> {"Time", "Count"}, 
      DefaultFont -> {"Helvetica-Bold", 12}, 
      Prolog -> AbsolutePointSize[4]]

Options[Regress]

report = 
      Regress[bactdata, {1, x}, x, 
        RegressionReport
           -> {SummaryReport, FitResiduals, SinglePredictionCITable, 
            ParameterConfidenceRegion}, ConfidenceLevel -> 0.90]

errors = FitResiduals /. report

ListPlot[errors, DefaultFont -> {"Helvetica-Bold", 12}]

Show[Graphics[ParameterConfidenceRegion /. report], Axes -> True, 
      AxesLabel -> {"Constant", "x"}, 
      DefaultFont -> {"Helvetica-Bold", 12}, AxesOrigin -> {180, -30}]

cookans = Regress[bactdata, {1, x}, x, RegressionReport -> {CookD}]

ListPlot[CookD /. cookans, DefaultFont -> {"Helvetica-Bold", 12}, 
      PlotRange -> {0, 1.60}, Prolog -> AbsolutePointSize[4], 
      AxesLabel -> {"Point", "CookD"}]

ListPlot[CookD /. cookans, DefaultFont -> {"Helvetica-Bold", 12}, 
      PlotRange -> {0, 0.3}, Prolog -> AbsolutePointSize[4], 
      AxesLabel -> {"Point", "CookD"}]

quadreport = 
      Regress[bactdata, {1, x, x^2}, x, 
        RegressionReport
           -> {SummaryReport, FitResiduals, SinglePredictionCITable, 
            ParameterConfidenceRegion}, ConfidenceLevel -> 0.90]

Probability

bintable2 = Table[Binomial[2, m], {m, 0, 8}];

bintable4 = Table[Binomial[4, m], {m, 0, 8}];

bintable6 = Table[Binomial[6, m], {m, 0, 8}];

bintable8 = Table[Binomial[8, m], {m, 0, 8}];

BarChart3D[{bintable8, bintable6, bintable4, bintable2}, 
      XSpacing -> 0.5, Boxed -> False, Axes -> False, 
      SolidBarEdgeStyle -> GrayLevel[0.5], 
      ViewPoint -> {2.94, -1.66, 0.245}]

geometric[p_] := Table[PDF[GeometricDistribution[p], x], {x, 0, 15}]

BarChart3D[{geometric[0.5], geometric[0.3], geometric[0.1]}, 
      Frame -> False, Boxed -> False, Ticks -> {None, Automatic, None}, 
      PlotLabel -> "Some Geometric Distributions", 
      AxesEdge -> {{1, -1}, {1, -1}, Automatic}, 
      AxesLabel -> {"p", "x", ""}, ViewPoint -> {3, 1, 1}]

Quantile[NormalDistribution[0, 1], 0.975]

N[Quantile[NormalDistribution[0, 1], 0.975], 22]

Quantile[NormalDistribution[354, 125], 0.935]

Quantile[StudentTDistribution[25.7], .975]

Quantile[StudentTDistribution[46], .975]

weibulldist = PDF[WeibullDistribution[α, β], x]

Mean[WeibullDistribution[α, β]]

Variance[WeibullDistribution[α, β]]

Plot[PDF[WeibullDistribution[4, 4], x], {x, 0, 8}, 
      AxesLabel -> {"X", "f"}, 
      DefaultFont -> {"Helvetica-Bold", 12}]

Plot[PDF[WeibullDistribution[1, 6], x], {x, 0, 8}, 
      AxesLabel -> {"X", "f"}, AxesOrigin -> {0, 0.04}, 
      DefaultFont -> {"Helvetica-Bold", 12}]

Plot[PDF[WeibullDistribution[2, 9], x], {x, 0, 18}, 
      AxesLabel -> {"X", "f"}, 
      DefaultFont -> {"Helvetica-Bold", 12}]

<< Graphics`Graphics3D`

Integrate[PDF[WeibullDistribution[2, 9], x], {x, 5, 14.6}]

weibullPDF[α_, x_] := 
      PDF[WeibullDistribution[α, 1], x]

weibullchart = 
      Plot3D[weibullPDF[α, x], {x, 10^(-1), 2}, {α, 
          0.5, 2.5}, PlotRange -> All, PlotPoints -> {30, 30}, 
        Ticks -> {Automatic, Automatic, None}, 
        PlotLabel -> "Weibull Distribution: density", 
        AxesEdge -> {Automatic, {1, -1}, Automatic}, 
        AxesLabel -> {"x", "  α", " "}]

cauchypdf[lambda_, x_] := PDF[CauchyDistribution[0, lambda], x]

cauchygraph = 
      Plot3D[cauchypdf[lambda, x], {x, -10, 10}, {lambda, 1, 4}, 
        PlotRange -> All, PlotPoints -> {40, 20}, 
        Ticks -> {Automatic, Automatic, None}, 
        PlotLabel -> "Some Cauchy Distributions", 
        AxesEdge -> {Automatic, {1, -1}, Automatic}, 
        AxesLabel -> {"x", "lambda", " "}]

ListPlot[Table[PDF[BinomialDistribution[12, .65], x], {x, 0, 12}], 
      AxesLabel -> {"X", "f"}, Prolog -> AbsolutePointSize[4], 
      DefaultFont -> {"Helvetica-Bold", 12}]

(∑_(x = 4)^8) PDF[BinomialDistribution[12, 0.65], x]

Bivariate Normal Distributions

k[r_] = 1/(2*Pi*Sqrt[1 - r^2])

f[x_, y_] := 
      k[r]*Exp[(-1/(2 (1 - r^2)))*(x^2 - 2*r*x*y + y^2)]

r = 1/2

Plot3D[f[x, y], {x, -3, 3}, {y, -3, 3}, 
      AxesLabel -> {"x", "y", "f"}, PlotPoints -> 40, 
      DefaultFont -> {"Helvetica-Bold", 12}]

ContourPlot[f[x, y], {x, -3, 3}, {y, -3, 3}, PlotPoints -> 40]

r = 0.7

Plot3D[f[x, y], {x, -3, 3}, {y, -3, 3}, 
      AxesLabel -> {"x", "y", "f"}, PlotPoints -> 40, 
      DefaultFont -> {"Helvetica-Bold", 12}, PlotRange -> {0, .25}]

ContourPlot[f[x, y], {x, -3, 3}, {y, -3, 3}, PlotPoints -> 40]

NIntegrate[f[x, y], {x, 0, 2}, {y, 0, 1}]

Race Car Problem

Race cars, numbered from 1 to 100, are running around a track. An observer notes the numbers on nine  cars. Assuming the cars appear in random order, what is the probability the median of the sample of nine cars is 35? We define a function giving the probability that the median of the sample is k.

f[k_] = Binomial[k - 1, 4]*Binomial[100 - k, 4]/Binomial[100, 9]

f[35]

% // N

ListPlot[Table[f[k], {k, 5, 95}], 
      DefaultFont -> {"Helvetica-Bold", 12}, 
      Prolog -> AbsolutePointSize[4]]

Generating Functions

We toss a fair die n times and examine the sums that occur as n varies.

g[t_] = (1/6)*Sum[t^n, {n, 1, 6}]

ListPlot[CoefficientList[(g[t]^2), t], 
      AxesLabel -> {"Sum", "f"}, Prolog -> AbsolutePointSize[4]]

ListPlot[CoefficientList[(g[t]^4), t], 
      AxesLabel -> {"Sum", "f"}, Prolog -> AbsolutePointSize[4]]

ListPlot[CoefficientList[(g[t]^12), t], 
      AxesLabel -> {"Sum", "f"}, Prolog -> AbsolutePointSize[4]]

Acceptance Sampling

From a lot of 1000 items, we draw a sample of size 100.  If the sample contains any defective items, we inspect the entire lot, replacing any defective items in the lot with good items. So if the lot contains d defective items before inspection, the lot shipped contains either d or 0 defective items. The average outgoing quality is the percentage of defective items sold and is given by

aoq[d_] = (d/1000)*Binomial[1000 - d, 100]/Binomial[1000, 100];

ListPlot[Table[aoq[d], {d, 0, 90}], AxesLabel -> {"D", "AOQ"}, 
      DefaultFont -> {"Helvetica-Bold", 12}, PlotRange -> {0, 0.0035}, 
      Prolog -> AbsolutePointSize[4]]

Recursions

We toss a coin, loaded to come up heads with probability p, until two heads occur in a row. 

Let a[n] denote the probability that the sequence HH occurs for the first time at the nth toss. Then

a[n_] := a[n] = (1 - p)*a[n - 1] + p*(1 - p)*a[n - 2]

a[1] = 0

a[0] = 0

a[2] = p^2

recurvals = Table[a[n], {n, 2, 10}] // Simplify

<< DiscreteMath`RSolve`

ans = 
      RSolve[{b[n] - (1 - p)*b[n - 1] - p*(1 - p)*b[n - 2] == 0, 
          b[0] == 0, b[1] == 0, b[2] == p^2}, b[n], n]

h[n_, p_] = b[n] /. ans[[1]]

h[2, 1/2]

-2 
      (((2 ((((1/4) - ((Sqrt[5])/4)))^2) (-(1/4) + ((Sqrt[5])/4)) 
              (-(1/2) + (1/2) (-(1/2) - ((Sqrt[5])/2))))/(Sqrt[5]))
               - ((2 (-(1/4) - ((Sqrt[5])/4)) ((((1/4) + ((Sqrt[5])/4)))^2) 
              (-(1/2) + (1/2) (-(1/2) + ((Sqrt[5])/2))))/(Sqrt[5])))

% // N

Table[h[n, 1/2], {n, 2, 10}] // Simplify

Table[h[n, 1/2]*2^n, {n, 2, 10}] // Simplify

\(\(














\)\)

German Tanks

samples = 
      Flatten[Table[{i, j, k}, {i, 1, 12}, {j, i + 1, 12}, {k, j + 1, 12}], 2]

estimates = 
      Table[(4/3)*Max[samples[[i]]] - 1, {i, 1, Length[samples]}]

sortest = Sort[estimates];

valuesest = Union[sortest]

frequest = 
      Table[Count[sortest, valuesest[[i]]], {i, 1, Length[valuesest]}]

Mean[estimates]

ListPlot[Transpose[{valuesest, frequest}], Frame -> True, 
      FrameLabel -> {"Estimate", "Frequency"}, 
      DefaultFont -> {"Helvetica-Bold", 12}, AxesOrigin -> {0, 0}]

And Finally

The object of computing is insight, not numbers.
		
Richard Hamming - Numerical Methods for Scientists and Engineers, Mc-Graw Hill Book Company, 1962