% eda12_13
% Atlantic Rocks dataset
% bootstrap error analysis for
% CaO/Na2O ratio of varimax factor 2
% note CaO/Na2O are elements 6/7

% load data
Draw = load('rocks.txt');
% 1 SiO2
% 2 TiO2
% 3 Al203
% 4 FeO-total
% 5 MgO
% 6 CaO
% 7 Na2O
% 8 K2O
Ns = size(Draw);
N = Ns(1);
M = Ns(2);

% number of realizations
Nr=1000;
ratio = zeros(Nr,1);

for ir = [1:Nr]

% randomply sample Draw
D = zeros(N,M);
D = Draw(unidrnd(N,N,1),:);

% compute factors and factor loadings using singular value decompostion
[U, S, V] = svd(D,0);

% keep only first five singular values
P=5;
F = V(:,1:P)';
C = U(:,1:P)*S(1:P,1:P);

% initial rotated factor matrix, FP, and rotation matrix, MR
MR=eye(P,P);
FP=F;

% spike these factors using the varimax procedure
k = [2, 3, 4, 5]';
Nk = length(k);

for iter = [1:3]
for ii = [1:Nk]
for jj = [ii+1:Nk]
    
% spike factors i and j
i=k(ii);
j=k(jj);

% copy factors from matrix to vectors
fA = FP(i,:)';
fB = FP(j,:)';

% standard varimax procedure to determine rotation angle q
u = fA.^2 - fB.^2;
v = 2* fA .* fB;

A = 2*M*u'*v;
B = sum(u)*sum(v);
top = A - B;

C = M*(u'*u-v'*v);
D = (sum(u)^2)-(sum(v)^2);
bot =  C - D;

q = 0.25 * atan2(top,bot);

cq = cos(q);
sq = sin(q);

fAp =  cq*fA + sq*fB;
fBp = -sq*fA + cq*fB;

% put back into factor matrix, FP
FP(i,:) = fAp';
FP(j,:) = fBp';

% accumulate rotation
mA = MR(i,:)';
mB = MR(j,:)';
mAp =  cq*mA + sq*mB;
mBp = -sq*mA + cq*mB;
MR(i,:) = mAp';
MR(j,:) = mBp';

end
end
end

% new factor loadings
CP=C*MR';

% ratio of elemnts 6 and 7 of factor 2
ratio(ir)=FP(2,6)/FP(2,7);

end

% compute histogram, convert to probability density distribution
Nb=100;
[rhist,bins] = hist(ratio,Nb);
Db=bins(2)-bins(1);
pbootstrap=rhist/(Db*sum(rhist));

% mean, variance
mr = Db * sum( bins .* pbootstrap );
vr = Db * sum( ((bins-mr).^2) .* pbootstrap );
disp('CaO/Na2O ratio of Factor 2');
disp(sprintf('mean %f variance %f',mr,vr));

% confidence intervals
Pbootstrap=Db*cumsum(pbootstrap);
ilo = find( Pbootstrap>=0.025, 1 );
ihi = find( Pbootstrap>=0.975, 1 );
rlo = bins(ilo);
rhi = bins(ihi);
disp(sprintf('%f < ratio < %f (95 percent confidence)',rlo,rhi));

figure(1);
clf;
set(gca,'LineWidth',2);
hold on;
axis( [bins(1), bins(Nb), 0, max(pbootstrap) ] );
plot( bins, pbootstrap, 'k-', 'LineWidth', 2);
plot( [mr, mr], [0, 0.05*max(pbootstrap)], 'k-', 'LineWidth', 2);
plot( [rhi, rhi], [0, 0.05*max(pbootstrap)], 'k-', 'LineWidth', 2);
plot( [rlo, rlo], [0, 0.05*max(pbootstrap)], 'k-', 'LineWidth', 2);
xlabel('CaO/Na2O ratio, r');
ylabel('p(r)');

