Principal Components - Monte Carlo

Overview


The principal components algorithm returns a set of orthonormal vectors along with a variance which represents the variance of that vector. (we assume that the data was centered prior to running the algorithm) To simulate the principal components, choose a distribution which represents the distribution of each component. (the most commonly used distribtuion is the normal distribtuion)

For each component you want to simulate, generate a random number from the distribution chosen for that component, with mean equal to zero and variance equal to the variance associated with that component.

Next multiply each component by the random number generated for that component and add the components together.

Definition


{% \vec{m} = \mathbb{E}(\vec{Z}) %}
{% cov(\vec{Z}) = \mathbb{E}[(A\vec{Z} - A\vec{m})(A\vec{Z} - A\vec{m})^T] %}
{% = \mathbb{E}[A(\vec{Z}-\vec{m})(\vec{Z}-\vec{m})^T A^T] %}
{% = A\mathbb{E}[(\vec{Z}-\vec{m})(\vec{Z}-\vec{m})^T]A^T %}
{% = A cov(\vec{Z}) A^T %}

Implementation





let la = await import('/lib/linear-algebra/v1.0.0/linear-algebra.mjs');
let mt = await import('/lib/statistics/moments/v1.0.0/moments.mjs');
let nm = await import('/lib/statistics/distributions/normal/v1.0.0/normal.mjs');
let data = [
	[1,2,3],
	[2,5,3],
	[6,3,4],
	[3,3,3],
	[1,4,2],
	[7,8,4],
];
let centered = mt.center(data);
let covariance = mt.covariance(centered);
let eig = la.eigenvectors(covariance);
let eigenvalues = eig.map(p=>p.eigenvalue);
let vectors = eig.map(p=>p.eigenvector);


The next step is to simulate a random vector. For this example, we will only use the first two principal components in the simulation. The distribution of each component is assumed to be normal. 


let nm = await import('/lib/statistics/distributions/normal/v1.0.0/normal.mjs');

//get the first vector
let v1 = vectors[0];
let v2 = vectors[1];

let rand1 = nm.random(0, eigenvalues[0]);
let rand2 = nm.random(0, eigenvalues[1]);

let result = v1.map((p,i)=>{
  return rand1*v1[i] + rand2*v2[i];
});
Try it!


Lastly we will compute a number of simulations and store them in a new array. Then compute the covariance of the new array of data. This new covariance matrix should look close to the covariance matrix of the original data. (Note, it will not be exact as we have ignored the third principal component) 


let data2 = [];
for(let i=0;i<500;i++){
  let rand1 = nm.random(0, eigenvalues[0]);
  let rand2 = nm.random(0, eigenvalues[1]);

  let result = v1.map((p,i)=>{
    return rand1*v1[i] + rand2*v2[i];
  });
  data2.push(result);
}

//for testing purposes
//this covariance matrix should look similar to the one computed above
let covar2 = mt.covariance(data2);
Try it!

Contents