Mathematical Principles for Scientific Computing and Visualization  
http://www.farinhansford.com/books/scv
by Gerald Farin and Dianne Hansford

This notebook demonstrates scatter  plots, correlations, line of regression, and dominant lines. These tools help us see trends in the data.

Input: two data sets. Here we have 8 weight and height pairs of data.
Output: mean, standard deviation, correlation coefficient, "time series" plots, scatter plot, the line of regression and the dominant line over the scatter plot. 

The dominant line is computed by revisiting PCA from Chapter 6.

Figures  12.1, 12.2, and 12.5 were created with this notebook.

See Chapter 12: Visualization of Empirical Data for details.

Updated May 2008

(* Load packages needed *)
Needs["PlotLegends`"]

(* Set output folder *)
SetDirectory[ToFileName[Extract["FileName"/.NotebookInformation[EvaluationNotebook[]],{1},FrontEnd`FileName]]];

(* Create data sets *)
weights={65,60,50,52,75,52,72,63};
heights={172,170,150,156,172,160,161,160};

(* Find mean of each data set *)
meanwt=N[Mean[weights]];
meanht=N[Mean[heights]];

(* Find standard deviation of data sets *)
stdwts=N[StandardDeviation[weights]];
stdhts=N[StandardDeviation[heights]];

(* Center each data set about the mean *)
mweights=weights-meanwt;
mheights=heights-meanht;

Print["Weights translated so mean is at the origin"]
Print[mweights]
Print[" "]
Print["Height translated so mean is at the origin"]
Print[mheights]
Print[" "]

(* Get the standard deviation of the modified data sets *)
(* And print it out *)
Print["Standard deviation of mweights = ",N[StandardDeviation[mweights]]]
Print["Standard deviation of mheights = ",N[StandardDeviation[mheights]]]
Print[" "]

(* Comput the correlation coefficient between the data sets *)
Print["Correllation coefficient for two data sets = ",Correlation[mweights,mheights]]
Print[" "]

(* FInd the slope of the line of regression *)
(* See equation (12.4) for details *) 
slope=((mweights.mheights)/(mweights.mweights));

Print["Line of regression: slope= ",slope]
Print[" "]


(* Set up for drawing the line of regression *)
linelst={{-15,-15 slope},{15,15 slope}};

(* Create a plot of the weights data set which has been translated so the it is centered about its mean *)
Print["Plot of mweights "]
picwts=ListPlot[mweights,AxesOrigin->{0,0},PlotStyle->PointSize[0.03`],PlotRange->{{0,9},{-15,15}},BaseStyle->{"FontFamily"->"Helvetica","FontSize"->11}]

(* Create a plot of the heights data set -- also modified to be centered about its mean *)
Print["plot of mheights "]
pichts=ListPlot[mheights,AxesOrigin->{0,0},PlotStyle->PointSize[0.03`],PlotRange->{{0,9},{-15,15}},BaseStyle->{"FontFamily"->"Helvetica","FontSize"->11}]

(* Create a table bringing together the weight and height data sets *)
(* This will be used for a scatter plot *)
joindata=Table[{mweights[[k]],mheights[[k]]},{k,1,Length[mweights]}];

Print["Weight and height data in one table as {weight,height} pairs: ",joindata]
Print[" "]

(* Create a scatter plot with the line of regression *)
Print["Scatter plot -- coordinates (weight,height) and add the line of regression"]
picscatter=ListPlot[{joindata, linelst},Joined->{False, True},PlotMarkers->Automatic,PlotStyle->{{PointSize[0.03]}, {Thickness[0.015]}}, BaseStyle->{"FontFamily"->"Helvetica","FontSize"->11},PlotRange->{{-15,15},{-15,15}}]

(* Output plots to postscript files *)
Export["timeseries_wts.eps",picwts];
Export["timeseries_hts.eps",pichts];
Export["scatterplot_wtshts.eps",picscatter];




(* Principal Components Analysis -- PCA *)

(* Copy the weight and height data *)
(* We have a table/matrix of (weight1,height1), (weight2, height2),.... *)
(* Having n rows and 2 columns *)
rmat=joindata
(* Create a transpose matrix having 2 rows and n columns *)
rtmat=Transpose[joindata]

(* Multiply the transpose matrix times the matrix *)
smat=Transpose[joindata].joindata;

Print["matrix from jointdata^T * joindata   smat = ",smat]
Print[" "]

(* Find the eigenvalues of the product matrix *)
Print["Eigenvalues of smat = ",Eigenvalues[smat]]
Print[" "]

(* Find the eigenvectors of the product matrix *)
(* The eigenvector corresponding to the largest eigenvalue is the dominant line of the data set *)
{e1,e2}=Eigenvectors[smat];

Print["Eigenvectors of smat:  e1=",e1,"  e2=",e2]
Print[" "]

(* Slope of the dominant line is the quotient of the components of the dominant eigenvector. Recall slope is rise/run  *)
domslope=((e1[[2]])/(e1[[1]]));

Print["Slope of dominant line = ",domslope]
Print[" "]

(* Create data to plot the dominant line *)
domline={{-15,-15 domslope},{15,15 domslope}};

(* Create a plot of the weight/height data, the dominant line, and the line of regression *)
Print["Plot of (weight,height) data, the dominant line, the line of regression"]
Print["Dominant line is brown. Line of regression is purple."]
picdominant=ListPlot[{joindata,linelst,domline},Joined->{False,True,True},PlotStyle->{{Blue,PointSize[0.03]},{Purple,Thickness[0.01`]},{Brown,Thickness[0.01`]}},PlotRange->{{-15,15},{-15,15}},BaseStyle->{"FontFamily"->"Helvetica","FontSize"->11}]



(* Output this plot to a postscript file *)
Export["scatterplot_dom.eps",picdominant]