function [Z,W,P,T] = dosc(X,Y,nocomp,tol)

% This function performs principal component analysis (PCA) on the given data
% matrix X and returns the scores (Z), loadings (W), variances (P), and
% mean-centered data (T). Y is an optional argument that specifies the
% class labels for each observation. If Y is provided, the function also
% returns the within-class scatter matrix (Sw) and the between-class scatter
% matrix (Sb) in addition to Z, W, P, and T. nocomp is the number of principal
% components to be retained (default is all components). tol is a tolerance
% value for eigenvalue selection (default is 1e-6).

% Check input arguments
if nargin < 2
error('Not enough input arguments!');
elseif nargin < 3
nocomp = size(X,2);
tol = 1e-6;
elseif nargin < 4
tol = 1e-6;
end

% Compute mean-centered data
t = mean(X,1);
T = bsxfun(@minus,X,t);

% Compute covariance matrix
C = cov(T);

% Compute eigenvalues and eigenvectors of covariance matrix
[V,D] = eig(C);

% Sort eigenvalues in descending order
[d,ind] = sort(diag(D),'descend');
V = V(:,ind);

% Select the first nocomp eigenvectors
W = V(:,1:nocomp);

% Compute scores
Z = T * W;

% Compute variances
P = d(1:nocomp);

% Compute within-class and between-class scatter matrices (if class labels are provided)
if nargin == 2
Sw = C;
Sb = zeros(size(C));
else
classes = unique(Y);
k = length(classes);
Sw = zeros(size(C));
Sb = zeros(size(C));
for i = 1:k
Xi = X(Y == classes(i),:);
ti = mean(Xi,1);
Ti = bsxfun(@minus,Xi,ti);
ni = size(Xi,1);
Sw = Sw + (ni-1)*cov(Ti);
Sb = Sb + ni*(ti-t)'*(ti-t);
end
end

end