For an assignment the Matlab code has to be rewritten into Python code.

formatcompact % For compact printing of results
mu = [23];
Cov = [21.5;1.54]
rng('default') % For reproducibility
R = mvnrnd(mu,Cov,100);


I rewrote it into:

import matplotlib.pyplot as plt
import math
from numpy.random import seed
from sklearn.preprocessing import scale # Data scaling
from sklearn.decomposition import PCA
import numpy as np
#mu = [2 3];
mu = np.array([2, 3])
#Cov = [2 1.5; 1.5 4];
Cov = np.array([[2, 1.5], [1.5, 4]])
#rng('default') % For reproducibility
#R = mvnrnd(mu,Cov,100);
R = np.random.multivariate_normal(mu, Cov, 100)


This line here

[loading_vector,score,latent,tsquared,variablity_explained,mu] = pca(R);


puzzles me. In my approach

R = scale(R.T)
pca = PCA(n_components=2, svd_solver='full')
pca.fit(R)

score = pca.transform(R)
eigvals, eigvect = np.linalg.eig(np.cov(R.T))


I am not familiar with the matlab counterpart (have not used it for some time), but I assume you are looking for the eigenvalues and eigenvectors of the covariance matrix from scikit-learn.

Since you already have the pca object and have fitted it to the data R, the values you are looking for are retrievable as object attributes:

loadings = pca.components_

This should return the same results as what you have by doing eigendecomposition on your covariance. The slight difference is that the scikit-learn pca counterpart would have already sorted both the squared and loadings by descending order of the singular value, whereas your np.linalg.eig would not have been sorted.