Statistical analysis with Mathematica - part 1

Load some standard statistics add-on packages

<<"LinearRegression`"

Example of use of a continuous distribution:

The probability density for value  x  in a gamma distribution is proportional to  (e^(-x/β))(x^(α-1))  for x>0,and is zero forx<0.

gdist = GammaDistribution[3,1]

Define the Cumulative Distribution Function for the above PDF :

cdfunction = CDF[gdist,x]

N[Evaluate[CDF[gdist,1]]]

Plot[cdfunction,{x,0,10}]

Obtain the PDF  by deriving the CDF :

gfunction = Dt[cdfunction,x]

Plot[gfunction,{x,0,10}]

Example of a simple linear fit to a set of data points:

points ={{0,-1},{1,2},{3,7},{4,6},{7,9},{11,14}};

plotpoints=ListPlot[points,PlotStyle->PointSize[0.02`]]

Fit the data points with a linear function a*(x + 2) :

Fit[points,x+2,x]

plotline=Plot[%,{x,-3,14}]

Show[plotpoints,plotline]

Let's perform the same fit with more details, using the Regress function of the Statistics`LinearRegression` add-on package:

regresult=Regress[points,{x+2},x,IncludeConstant->False, RegressionReport->{FitResiduals,SinglePredictionCITable,ParameterConfidenceRegion}]

Another choice for the Regression Report :

regmore = Regress[points,{x+2},x,IncludeConstant->False,RegressionReport->{ParameterTable}]

errors = FitResiduals /. regresult

Use the Transpose built-in function to extract single columns from the Single Prediction Confidence Interval table:

{observed,predicted,se,ci}=Transpose[(SinglePredictionCITable /. regresult)[[1]]]

Pair predicted, lower and upper confidence limits with the corresponding x values:

xval = Map[First,points];predicted = Transpose[{xval,predicted}]; lowerCI=Transpose[{xval,Map[First,ci]}]; upperCI=Transpose[{xval,Map[Last,ci]}];
errors=Transpose[{xval,errors}];

Display data, fitted curve, 95% confidence interval for predicted responses of single observations:

ListPlot[{points,predicted,lowerCI,upperCI},PlotMarkers->Automatic,PlotStyle->{Automatic,Automatic,Dashing[{0.02`,0.02`}],Dashing[{0.02`,0.02`}]},Joined->{False,True,True,True}]

Display fit residuals:

ListPlot[{errors},Filling->Axis]

Let's now see an example of a parabolic fit (only constant term + quadratic term).

data = {{0.05,89},{0.097,98},{0.124,109},{0.175,144},{0.201,156},{0.247,180},{0.299,209},{0.347,251}};

ListPlot[{data},PlotMarkers->Automatic]

regres2 = Regress[data, {1,x^2},x,RegressionReport->{FitResiduals,SinglePredictionCITable,ParameterTable,ParameterConfidenceRegion}]

{params2,stderr2,tstat2,pval2}=Transpose[(ParameterTable /. regres2)[[1]]]

errors2 = FitResiduals /. regres2

{obs2,pred2,se2,ci2}=Transpose[(SinglePredictionCITable /. regres2)[[1]]];

xval2 = Map[First,data];pred2 = Transpose[{xval2,pred2}]; lowerCI2=Transpose[{xval2,Map[First,ci2]}]; upperCI2=Transpose[{xval2,Map[Last,ci2]}];
errors2=Transpose[{xval2,errors2}];

Display data,fitted curve,95% confidence interval for predicted responses of single observations:

ListPlot[{data,pred2,lowerCI2,upperCI2},PlotMarkers->Automatic,PlotStyle->{Automatic,Automatic,Dashing[{0.02`,0.02`}],Dashing[{0.02`,0.02`}]},Joined->{False,True,True,True}]

Display fit residuals:

ListPlot[{errors2},Filling->Axis]

Display the parameters' best values and the 95% confidence region for the two parameters:

Show[Graphics[ParameterConfidenceRegion /. regres2],Graphics[Point[params2]],AspectRatio->1,Axes->True, AxesLabel->{"constant","x squared"}]

Finally, let' s perform the fit with a full 2nd degree polynomial :

regres3 = Regress[data, {1,x,x^2},x,RegressionReport->{FitResiduals,SinglePredictionCITable,ParameterTable,ParameterConfidenceRegion}]