close all
clear
clc
formatSpec = '%f';
sizeA = [1 Inf];
x=1:1001;

M  = 2;

fileID1 = fopen('first_process.txt','r');
realization_mat = fscanf(fileID1,formatSpec,sizeA);
realization_mat=reshape(realization_mat, 50, 1001);
realization_mat=realization_mat';
IndividualObservations=size(realization_mat, 1);

fileID2 = fopen('second_process.txt','r');
realization_mat2 = fscanf(fileID2,formatSpec,sizeA);
realization_mat2=reshape(realization_mat2, 50, 1001);
realization_mat2=realization_mat2';


realization_mat_merged=cat(1, realization_mat-mean(realization_mat, 2), realization_mat2-mean(realization_mat2, 2));
realization_mat_merged=realization_mat_merged';

[coeff, score, latent, tsquared, explained, mu]=pca(realization_mat_merged', 'Centered',false);

reconstruccion_aPC=score*coeff'+ repmat(mu,2002,1);
reconstruccion1_aPC=reconstruccion_aPC(1:1001, :)+mean(realization_mat, 2);
reconstruccion2_aPC=reconstruccion_aPC(1002:2002, :)+mean(realization_mat2, 2);


N1 = IndividualObservations;
w = [.005;ones(999,1)*.001;.005];

C=cov(realization_mat_merged);

    
    
[U, d] = eig(C);        
d = diag( d ); 

U= real(U);
d = real(d); 

d = d(end:-1:1); 
U = U(:,end:-1:1);


rDim = 50;
randVar=U'*realization_mat_merged';

reconstruction=randVar'*U';
reconstruction=reconstruction';

reconstruction1=reconstruction(1:N1, :);
reconstruction2=reconstruction(1+N1:end, :);

reconstruction1=reconstruction1+mean(realization_mat, 2);
reconstruction2=reconstruction2+mean(realization_mat2, 2);

figure
subplot(2,1,1)
plot(score(:,1:4))
title('scores')
subplot(2,1,2)
plot((sqrt(d(1))*U(:,1))*(-1))
hold on
for i=2:4
    plot(sqrt(d(i))*U(:,i))
end
title('evec')


figure
subplot(2,2,1)
histogram(randVar(1,:))
title('RV 1 handmade')
subplot(2,2,2)
histogram(randVar(2,:))
title('RV 2 handmade')
subplot(2,2,3)
histogram(randVar(3,:))
title('RV 3 handmade')
subplot(2,2,4)
histogram(randVar(4,:))
title('RV 4 handmade')


figure
subplot(2,2,1)
histogram(coeff(:,1))
title('RV 1 coeff')
subplot(2,2,2)
histogram(coeff(:,2))
title('RV 2 coeff')
subplot(2,2,3)
histogram(coeff(:,3))
title('RV 3 coeff')
subplot(2,2,4)
histogram(coeff(:,4))
title('RV 4 coeff')



