function [latVars, PCs, What, Up, Dp] = NOlatentVariables(numLatVar, numObsVar, sampleSize, upperVar)
%NOlatentVariables A simulation program to show if a PCA can capture non-orthogonal latent variables.
%
% USE:
% [latVars, PCs, What, Up, Dp] = NOlatentVariables(numLatVar, numObsVar, sampleSize)
%
% INPUTS:
% numLatVar = the number of latent variables
% numObsVar = the number of observable variables
% sampleSize = the number of samples
% upperVar = a numLatVar x 1 vector.  These values define the upper limit of the variance
%            for each latent variable.  
% 
% OUTPUTS:
% latVars = A matrix of latent variables (sampleSize x numLatVar)
% PCs = A matrix of principal component scores (sampleSize x numObsVar)
% What = Observables that were created from latent variables. 
%
% Steps:
% 1. Create an arbitrary Up, which contains non-orthogonal vectors.
% 2. Create latent variables (Y).
% 3. Create a data set (What) from these latent variables What = Y x Up'.
% 4. Add an error term to What: W2hats = What + error.
% 5. Conduct PCA on W2hats.
% 6. Compare PCA scores to the latent variables. 

% Tomo Eguchi
% 22 February 2002 

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

% create variance matrix:
Dp = zeros(numLatVar);
for i = 1:numLatVar,
   Dp(i, i) = random('unif', 1, upperVar(i), 1, 1);
end

% create arbitrary mixing matrix:
Up = random('norm', 3, 1, numObsVar, numLatVar);

% create latent variables:
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