11
$\begingroup$

There is a larger matrix (1500 rows x 40 columns), 1500 observations x 40 variables. then I follow the procedures of PCA(Principal components analysis), 1. find correlation 2. find eigenvalues 3. find Eigenvectors I have got the result (Mathematica)

dataCorrEigenvalues/Total@dataCorrEigenvalues

{0.647833, 0.128731, 0.0843738, 0.0519215, 0.0246577, 0.018331, \
0.0100494, 0.00657219, 0.0054721, 0.00373078, 0.00310175, 0.00244999, \
0.0022292, 0.00190861, 0.00166728, 0.00124446, 0.00113064, \
0.00093684, 0.000673087, 0.000579798, 0.00049716, 0.000425554, \
0.000371012, 0.000261027, 0.000225517, 0.000173631, 0.000133479, \
0.000128954, 0.000103792, 0.0000853669}

FoldList[Plus, dataCorrEigenvalues/Total@dataCorrEigenvalues]

{0.647833, 0.776564, 0.860938, 0.91286, 0.937517, 0.955848, 0.965898, \
0.97247, 0.977942, 0.981673, 0.984775, 0.987225, 0.989454, 0.991362, \
0.99303, 0.994274, 0.995405, 0.996342, 0.997015, 0.997595, 0.998092, \
0.998517, 0.998888, 0.999149, 0.999375, 0.999548, 0.999682, 0.999811, \
0.999915, 1.}

As I know the first 6 components explain about 95% of the variability, However, I don't understand how to use these components for data analysis.

The projection Matrix "w" as the following

w = dataCorrEigenvectors[[All, 1 ;; 4]];

I try to calculate the dot product of W against data, however, it seems the result is not same as expectation.

PC5 = data.w;

Please feel free to command and advise what I should do. Thank you.

$\endgroup$
  • 2
    $\begingroup$ I think this question about Principal Components is better suited for Cross Validated where questions about statistics are addressed as I don't see a specific Mathematica question. If the question is transferred to that site, I would add a sentence or two about your objectives. For example, do you plan to perform a regression? Are you looking to interpret the principal components? $\endgroup$ – JimB Jan 21 '16 at 15:08
  • 1
    $\begingroup$ Transferred to Cross Validated stats.stackexchange.com/questions/191808/… $\endgroup$ – madeinQuant Jan 21 '16 at 15:45
  • 1
    $\begingroup$ Nelson, you should also realize that PCA is built-in to Mathematica with the PrincipalComponents function. I do agree with you though that it's highly educational to be able to do it by hand. I think the only problem in your approach might be that you are not standardizing your data nor are you transposing the eigenvector matrix in accordance with how Mathematica defines an eigenvector problem. In other words, try Standardize[data].Transpose[dataCorrEigenvectors] instead and see how you fare. $\endgroup$ – MarcoB Jan 21 '16 at 17:12
  • $\begingroup$ Dear MacoB, I would try it again. Thank you. $\endgroup$ – madeinQuant Jan 22 '16 at 0:41
21
$\begingroup$

I generated a multivariate data set to play with called fakedata (see at the bottom of the post). This contains five groups of samples, each in 25-replicate. Each sample is described by five descriptors with very different ranges (see at the bottom of the post for the generation of the data). fakedata therefore is a $125\times5$ matrix: this is smaller than yours, but these methods are obviously scalable.

First of all, as I mentioned in the comments, you should be aware of the fact that Mathematica has a built-in PrincipalComponents function

builtinscores = PrincipalComponents[fakedata, Method -> "Correlation"];

Having said that, however, trying to roll your own is simple and instructive. In addition, to the fun of doing it by hand, it also gives you access to the data you would need to construct e.g. loading plots and other interesting descriptors of the multivariate analysis.

As you mentioned, one obtains the transformation matrix by calculating the eigenvectors of the correlation matrix of the original data. In order to obtain the transformed PCA scores from the original data, you need to Standardize your data first, then dot-multiply it by those eigenvectors (check out Mathematica's definition of Eigensystem to see why we want to Transpose those eigenvectors first):

manscores = Standardize[fakedata].Transpose[Eigenvectors[Correlation[fakedata]]];

... and this is really all there is to it!


In fact, we can graphically check that we obtain the same PCA scores with the two methods. Here is a quick helper function to plot a 2D PCA score plot:

Clear[plotscores]
plotscores[scores_] := ListPlot[
   Partition[scores[[All, 1 ;; 2]], 25],
   PlotRange -> All, PlotRangePadding -> Scaled[0.1],
   PlotMarkers -> {Automatic, Medium}, PlotStyle -> ColorData[22, "ColorList"],
   AspectRatio -> 1, Frame -> True, Axes -> False, ImageSize -> Scaled[0.33]
 ]

Using this function, we can generate score plots from the scores calculated before by hand vs using the built-in function:

builtin = plotscores[builtinscores];
manual = plotscores[manscores];
Grid[{Style[#, FontSize -> 18] & /@ {"Built-in", "Manual"}, {builtin, manual}}]

manual vs. built-in

You will notice that the two plots are identical up to a sign change. Of course sign is irrelevant here, so we can adjust the sign of the ordinate of the manually calculated:

manualadj = plotscores[Cases[manscores, {x_, y_, rest__} -> {x, -y, rest}]];
Grid[{Style[#, FontSize -> 18] & /@ {"Built-in", "Manual (adj)"}, {builtin, manualadj}}]

sign-changed manual scores vs. built-in


As a small aside, in order to calculate percentage contributions of each principal component you can simply use:

Normalize[Eigenvalues[Correlation[fakedata]], Total]

(* Out: {0.456492, 0.261667, 0.202659, 0.0648009, 0.0143812} *)

For cumulative contributions, Accumulate is your friend, rather than using your FoldList expression:

Accumulate[%]

(* Out: {0.456492, 0.718159, 0.920818, 0.985619, 1.} *)

Finally, here is how I generated fakedata:

Clear[samplereplicates]
samplereplicates[mu_ /; (Dimensions[mu] == {5}), n_Integer] :=
 RandomVariate[#, n] &@
  MultinormalDistribution[
   mu,
   {{0.0001, 0., 0., 0., 0.},
    {0., 0.0004, 0., 0., 0.},
    {0., 0., 2250000, 0, 0.},
    {0., 0., 0, 5290000, 0.},
    {0., 0., 0., 0., 0.000025}}
  ]

SeedRandom[20160122]
fakedata = Join[
   samplereplicates[{0.21, 0.85, 28000, 15000, 0.185}, 25],
   samplereplicates[{0.15, 0.65, 35000, 75000, 0.185}, 25],
   samplereplicates[{0.35, 0.80, 25000, 25000, 0.185}, 25], 
   samplereplicates[{0.02, 0.50, 45000, 50000, 0.185}, 25],
   samplereplicates[{0.46, 0.70, 55000, 15000, 0.185}, 25]
 ];
$\endgroup$
  • $\begingroup$ Dear MarcoB, Thanks for you advice. I was using the python package scikit-learn for machine learning. However, I hope I can figure out the internal mechanics and wants to implement it by hand to understand it better. $\endgroup$ – madeinQuant Jan 25 '16 at 4:30
  • $\begingroup$ @NelsonMok You are welcome. Hope this helps! $\endgroup$ – MarcoB Jan 25 '16 at 4:32
  • 1
    $\begingroup$ @NelsonMok If you think that my answer fulfills your requirements, I would appreciate it if you would consider accepting it officially (you can click on the grey check mark next to the answer to do so). $\endgroup$ – MarcoB Jan 26 '16 at 1:40
  • $\begingroup$ I would construct a PCA index; the eigenvectors of fakedata {0.456492, 0.718159, 0.920818, 0.985619, 1.} The first component already explains about 46% of the variability in the 125 observations with 5 variables. So, I did a dot product (fakedata x first-component) and I get a PCA index (125x1). However, 46% is not enough, at least 93% of the variability. How to construct a PCA index with more component and to make a PCA index is accuracy? Please feel free to comment. $\endgroup$ – madeinQuant Jan 27 '16 at 8:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.