% eda08_07
%
clear all;

% Part 1: Make synthetic chemical data with a wide range of concentrations

% typical analytic errors
%    [Fe        Cu        Zn        Pb        As        Ag        Au]
se = [0.010000, 0.010000, 0.010000, 0.001000, 0.001000, 0.0000010, 0.0000001 ]';

% 7 elements, three factors
M = 7;
P = 3;
%    [Fe        Cu        Zn        Pb        As        Ag        Au]
v1 = [0.200000, 0.004000, 0.010000, 0.010000, 0.001000, 0.0000100, 0.0000010 ]';
v2 = [0.050000, 0.180000, 0.180000, 0.030000, 0.000200, 0.0000400, 0.0000005 ]';
v3 = [0.100000, 0.080000, 0.110000, 0.010000, 0.000700, 0.0000100, 0.0000001 ]';

F = [v1, v2, v3]';

% randomly chosen loadings
N = 10000;
c1 = random('uniform',0,1,N,1);
c2 = random('uniform',0,1,N,1);
c3 = 1-(c1+c2);
C = [c1, c2, c3];

% true sample matrix
Strue = C*F;

% observed sample matrix is true plus Normally-distributed random noise
for i=[1:M]
    Sobs(:,i) = Strue(:,i) + random('Normal',0,se(i),N,1);
end

Nbins = 100;
Fe_bins = 5*[0:Nbins-1]'/(Nbins-1);
Au_bins = 5*[0:Nbins-1]'/(Nbins-1);

% Part two factor analysis with no weighting

[U1,SIGMA1,V1] = svd(Sobs,0);
Fpre1 = V1'; % factors
Cpre1 = U1*SIGMA1;    % loadings

% P=2 significant factors
Fpre1P = Fpre1(1:P,:);
Cpre1P = Cpre1(:,1:P);
Spre1P = Cpre1P*Fpre1P;

figure(1);
clf;

% error
e1 = abs((Strue-Spre1P));
Fe1hist = hist(e1(:,1)./se(1),Fe_bins);
Au1hist = hist(e1(:,M)./se(M),Au_bins);
subplot(1,2,1);
hold on;
set(gca,'LineWidth',2);
axis( [Fe_bins(1), Fe_bins(end), 0, 1] );
plot( Fe_bins, Fe1hist/max(Fe1hist), 'k-', 'LineWidth', 2 );
xlabel('error/se');
ylabel('count');
subplot(1,2,2);
set(gca,'LineWidth',2);
hold on;
axis( [Au_bins(1), Au_bins(end), 0, 1] );
plot( Au_bins, Au1hist/max(Au1hist), 'k-', 'LineWidth', 2 );
xlabel('error/se');
ylabel('count');

% Part three factor analysis with weighting by analytic precision
w = 1./se;
[U2,SIGMA2,V2] = svd(Sobs*diag(w),0);
Fpre2 = V2'*diag(1./w); % factors
Cpre2 = U2*SIGMA2;    % loadings

% factors
Fpre2P = Fpre2(1:P,:);
Cpre2P = Cpre2(:,1:P);
Spre2P = Cpre2P*Fpre2P;

% error
e2 = abs((Strue-Spre2P));
Fe2hist = hist(e2(:,1)./se(1),Fe_bins);
Au2hist = hist(e2(:,M)./se(M),Au_bins);
subplot(1,2,1);
plot( Fe_bins, Fe2hist/max(Fe2hist), '-', 'LineWidth', 2,  'Color', [0.8,0.8,0.8] );
subplot(1,2,2);
plot( Au_bins, Au2hist/max(Au2hist), '-', 'LineWidth', 2, 'Color', [0.8,0.8,0.8] );


