function [T, P, X_res, X_orig] = pca_nipals(X,numpcs)

if argn(2) < 2 then
    numpcs = rank(X);
end

if numpcs > rank(X) then
    numpcs = rank(X);
end

format(16);

X_orig=X;

[sam,var]=size(X);

P=zeros(var,numpcs);//memory allocation (Scilab) for the loadings matrix
T=zeros(sam,numpcs);//memory allocation (Scilab) for the scores matrix
X_res=zeros(sam,var);//memory allocation (Scilab) for the residual matrix

threshold=1e-6;// threshold constant for convergence test

for j=1:numpcs

t=X(:,j);
t_pre=t-t;
        
        while (norm(t-t_pre)) > threshold
        p=(X'*t)/(norm(X'*t));//calculate and normalise the loading vector
        t_pre=t;//saving the original scores vector into  t_pre
        t=X*p;//calculating the new score vector (product of eigenvalue and eigenvector)
        end

X=X-t*p';//X- will be equal to the residual matrix (Error), X_orig= original input matrix
T(:,j)=t;//extracting score vectors
P(:,j)=p;//extracting loading vectors

end

// Eigenvalue and Variance Calculations

Eig=diag(T'*T); //Eigenvalues
Total_eig=sum(Eig);

for k=1:numpcs

    Eig_val_perc(k) =(Eig(k)/Total_eig)*100; // Calculating Variation percentage per component

end

endfunction



