function [X, Y, PCx, U, eigVord, lambda] = findModel(n, numVar)
% findModel A simulation program for recreating data set from a covariance matrix.
%
% This program recreate a data set by using a covariance matrix.
% Algorithm:
% 1. Define a covariance matrix (covX)
% 2. Conduct a spectral decomposition of the covariance matrix. covX = UDU'
% 3. Create a data set Y with means 0 and variances according to the covarianc matrix
%    It uses the eigen values of covX.
% 4. Create an observable data set X = Y * U'.
% 5. Conduct PCA on X.
% 6. Plot PC versus Y, which should be 1-to-1 plots.
%
% USE:
% [X, Y, PCx, U, eigVord, lambda] = findModel(n, numVar)
%
% INPUT:
% n = the number of observations.
% numVar = the number of variables.
%
% OUTPUT:
% X = data matrix. X = Y * U'.
% Y = data set with mean zero and variance = D.
% PCx = principal component scores.
% U = a matrix of eigen vectors of covariance matrix of X.
% eigVord = a matrix of eigenvectors of X.
% lambda = a diagonal matrix of eigenvalues of covariance matrix of X.
%
% Note that X is created from the covariance matrix.  

% Tomo Eguchi
% 13 February 2002

% create the covariance matrix of X.
% first create psuedo data from normal distributions.
% this is only for creating the covariance matrix of X
pData = zeros(n, numVar);
for i = 1:numVar,
   mean = unif(100, 200);   % mean
   sd   = unif(10, 100);   % sd
   pData(:, i) = normrnd(mean, sd, n, 1);  % create data.
end
SigmaX = cov(pData);

% Analyses start here *************************************************** 
% spectral decomposition of SigmaX and order eigenvalues and eigenvectors:
[U, lambda] = eigOrder(SigmaX);
Y = zeros(n, numVar);          % initialize the Y matrix

% Y contains Gaussian random variables with means 0 and variances = eigenvalues. 
for i = 1:numVar,
   Y(:, i) = normrnd(0, sqrt(lambda(i,i)), n, 1);
end

% create X, i.e., observables according to the covariance matrix of X.
X = Y * U';

% conduct PCA on X:
[PCx, eigVord, eigDord] = pca(X);

% End of analyses *******************************************************

% Printout results
fprintf(1, '\n Sample size of %d\n', n);
fprintf(1, '\n Covariance matrix of X was:\n')
for i = 1:numVar,
   for j = 1:numVar,
      fprintf(1, '   %.2f', SigmaX(i, j));
   end
   fprintf(1, '\n');
end

% Statistics:
covX = cov(demean(X));
fprintf(1, '\n Covariance Matrix of observed X (demeaned) was:\n')
for i = 1:numVar,
   for j = 1:numVar,
      fprintf(1, '   %.2f', covX(i, j));
   end
   fprintf(1, '\n');
end

fprintf(1, '\n Eigenvectors of the original cov(X) were:\n');
for i = 1:numVar,
   for j = 1:numVar,
      fprintf(1, '   %.4f', U(i, j));
   end
   fprintf(1, '\n');   
end

fprintf(1, '\n Eigenvectors of the sample cov(X) were:\n');
for i = 1:numVar,
   for j = 1:numVar,
      fprintf(1, '   %.4f', eigVord(i, j));
   end
   fprintf(1, '\n');   
end

fprintf(1, '\n Eigenvalues of the original cov(X) were:\n');
for i = 1:numVar,
   fprintf(1, '   %.4f', lambda(i,i));
end
fprintf(1, '\n');

fprintf(1, '\n Eigenvalues of the sample cov(X) were:\n');
for i = 1:numVar,
   fprintf(1, '   %.4f', eigDord(i,i));
end
fprintf(1, '\n');

% plot Y and PC scores to see if they match
for i = 1:numVar,
   plot(Y(:, i), PCx(:,i), '.');
   Xlabel = ['Y', num2str(i)];
   Ylabel = ['PC', num2str(i)];
   xlabel(Xlabel);
   ylabel(Ylabel);
   hold on;
   pause;
   hold off;
end