function [latVars, PCs] = latentVariables(numLatVar, numObsVar, sampleSize)
%latentVariables A simulation program to show if a PCA can capture latent variables.
%
% USE:
% [latVars, PCs] = latentVariables(numLatVar, numObsVar, sampleSize)
%
% INPUTS:
% numLatVar = the number of latent variables
% numObsVar = the number of observable variables
% sampleSize = the number of samples
% 
% OUTPUTS:
% latVars = A matrix of latent variables (sampleSize x numLatVar)
% PCs = A matrix of principal component scores (sampleSize x numObsVar)
%
% Steps:
% 1. Create an arbitrary covariance matrix (SigmaW: numObsVar x numObsVar)
% 2. Spectral decomposition of SigmaW = UDU'
% 3. Create the partition of U (Up), by using the first two eigen values and eigen vectors.
% 4. Create latent variables (Y) by using a couple of eigen values.
% 5. Create a data set (What) from these latent variables What = Y x Up'.
% 6. Add an error term to What: W2hats = What + error.
% 5. Conduct PCA on W2hats.
% 6. Compare PCA scores to the latent variables. 

% Tomo Eguchi
% 18 February 2002 started
% 19 February 2002 

%numLatVar = 2; numObsVar = 5; sampleSize = 5000;

% create covariance matrix:
covW = createCov(numObsVar, random('unif',10, 30, 1, 5));

% spectral decomposition of covW
[U, D] = eigOrder(covW);

% pick up first few.
Up = U * eye(numObsVar, numLatVar);
Dp = eye(numLatVar, numObsVar) * D * eye(numObsVar, numLatVar);

% create latent variables from eigen values:
latVars = cov2data(Dp, sampleSize);

% Create observables:
What = latVars * Up';

% add error term
W2hats = What + random('norm', 0, 1, sampleSize, numObsVar);  

% Conduct PCA on W2hats
[PCs, eigenVord, eigenDord] = pca(W2hats);

% Plot results:
for k = 1:numLatVar,
   figure(k)
   plot(PCs(:, k), latVars(:, k), '.');
   xstr = ['PC', num2str(k)];
   ystr = ['Latent Variable ', num2str(k)];
   xlabel(xstr);
   ylabel(ystr);
end