function [train_pca,test_pca,dataset_cumsum,percent_explained] = pcaForRF(train,test,threshold)

% This function performs PCA on the training dataset and applies the same
% transformation to the testing dataset. It returns the transformed
% datasets, cumulative sum of variance explained by each principal
% component, and the percentage of variance explained by each principal
% component.
%
% Inputs:
% train - Training dataset with observations in rows and features in
% columns.
% test - Testing dataset with observations in rows and features in columns.
% The number of columns must match the number of columns in the
% training dataset.
% threshold - A threshold value (between 0 and 1) that determines the
% number of principal components to keep. The function will
% keep the minimum number of principal components required
% to explain the threshold fraction of the variance in the
% dataset.
%
% Outputs:
% train_pca - Transformed training dataset.
% test_pca - Transformed testing dataset.
% dataset_cumsum - Cumulative sum of variance explained by each principal
% component.
% percent_explained - Percentage of variance explained by each principal
% component.

% Compute mean and standard deviation of training data
train_mean = mean(train);
train_std = std(train);

% Standardize the training and testing data
train_stdz = (train - train_mean) ./ train_std;
test_stdz = (test - train_mean) ./ train_std;

% Compute covariance matrix of the standardized training data
cov_matrix = cov(train_stdz);

% Compute eigenvectors and eigenvalues of the covariance matrix
[eig_vectors, eig_values] = eig(cov_matrix);

% Sort the eigenvectors in descending order of eigenvalues
[eig_values, idx] = sort(diag(eig_values), 'descend');
eig_vectors = eig_vectors(:, idx);

% Compute cumulative sum of variance explained by each principal component
variance_explained = eig_values / sum(eig_values);
dataset_cumsum = cumsum(variance_explained);

% Compute number of principal components required to explain the threshold
% fraction of the variance in the dataset
num_components = find(dataset_cumsum >= threshold, 1, 'first');

% Compute percentage of variance explained by each principal component
percent_explained = variance_explained * 100;

% Transform the standardized training and testing data using the
% eigenvectors
train_pca = train_stdz * eig_vectors(:, 1:num_components);
test_pca = test_stdz * eig_vectors(:, 1:num_components);