N=90;           %%% Number of subjects
kys=6;          %%% Number of items per factor, should be even
factN=5;        %%% Number of factors, should be 5
minusOne = 1;   %%% =1 if there are no inversions

fStr=sqrt(0.6); %%% The number is the s2n ratio

f1str = fStr;
f2str = fStr;
f3str = fStr;
f4str = fStr;
f5str = fStr;
%%% Factor strengths

Bmaara=kys/2; 

target  = cat(1, ones(Bmaara,1)*temp(1,:), minusOne*ones(Bmaara,1)*temp(1,:), ...
ones(Bmaara,1)*temp(2,:), minusOne*ones(Bmaara,1)*temp(2,:),  ...
ones(Bmaara,1)*temp(3,:), minusOne*ones(Bmaara,1)*temp(3,:),  ...
ones(Bmaara,1)*temp(4,:), minusOne*ones(Bmaara,1)*temp(4,:),  ...
ones(Bmaara,1)*temp(5,:), minusOne*ones(Bmaara,1)*temp(5,:));
%%% Rotation target matrix


AcorKys  = cat(1, f1str*ones(Bmaara,1)*temp(1,:), f1str*minusOne*ones(Bmaara,1)*temp(1,:), ...
f2str*ones(Bmaara,1)*temp(2,:), f2str*minusOne*ones(Bmaara,1)*temp(2,:), ...
f3str*ones(Bmaara,1)*temp(3,:), f3str*minusOne*ones(Bmaara,1)*temp(3,:), ...
f4str*ones(Bmaara,1)*temp(4,:), f4str*minusOne*ones(Bmaara,1)*temp(4,:), ...
f5str*ones(Bmaara,1)*temp(5,:), f5str*minusOne*ones(Bmaara,1)*temp(5,:));
%%% Signal variance

Ekys  = ones(N, kys*factN);  %%% Error variance

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

mainCounter=[0 0 0 0 0];  
selitys=0; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for mainloop=1:20
    
%%%%%%%%%%%%  Data generation
    
    F=random('Normal', 0, 1, [N factN]);   %%% Factor scores
    E=random('Normal', 0, Ekys);  %%% Error
    
    K = F*AcorKys' + E;   %%% Generate K

    [L2,scores,latent]=princomp(K); %% PCA; K(columns centered) = scores*L'
    R = L2(:,[1:factN]);  %%% Retain 5 first components

%%%%%%%%%%%%  Analyysi

    [d1,Z1,tr1] = procrustes(target, R);   %%% Procrustes rotation

    mainCounter(mainloop,:)=dot(target,Z1)./sqrt(dot(target,target).*dot(Z1,Z1)); %%% Save congruence 
    selitys(mainloop)=sum(latent([1:factN]))/sum(latent);  %%% Save variance explained
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

mainCounter=cat(1,mainCounter(:,1), mainCounter(:,2), mainCounter(:,3), mainCounter(:,4), mainCounter(:,5));
%%%  Pooling data over factors

[num2str(fStr^2), ' & ', num2str(mean(selitys),3), ' & ',  num2str(prctile(mainCounter,50),3), ' & ', num2str(prctile(mainCounter,95),3), ...
    ' & ',num2str(prctile(mainCounter,99),3)]
%%%  Print results