lab_pca.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c75558614ccda6a5",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "### Inverse Problems and Parameter Estimation, GEOS 627/427, University of Alaska Fairbanks\n",
    "\n",
    "- script lab_pca.ipynb\n",
    "- the goal of this tutorial is to establish/review the connections among eigen decomposition,\n",
    "singular value decomposition, principal component analysis, covariance matrix, correlation matrix, etc.\n",
    "\n",
    "- two key operations to the columns of the data matrix\n",
    "(1) centered: subtract mean\n",
    "(2) standardized (or scaled): divide by the standard deviation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cfae8e9b51e13d8b",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c66242ba-ed62-48e9-bdcb-fbca8583562e",
   "metadata": {},
   "outputs": [],
   "source": [
    "%run lib_header.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f635e43f-b512-499d-b726-22482967fb24",
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    import scikit\n",
    "except ImportError:\n",
    "    import os\n",
    "    os.system('pip install scikit-learn')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d02f906",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append('data/')\n",
    "import pandas as pd\n",
    "from pandas.plotting import scatter_matrix\n",
    "from scipy import linalg as LA\n",
    "from scipy import stats\n",
    "from scipy.stats import zscore\n",
    "from sklearn.decomposition import PCA\n",
    "\n",
    "from lib_inverse import *                 # svdmat, pearsonr_ci\n",
    "from load_pca_data import load_pca_data   # this resides in the ./data/ folder\n",
    "from lib_geos import corrcov"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c65f1991",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load data\n",
    "# idata = 1: City ratings data (Matlab tutorial)\n",
    "# idata = 2: Protein consumption data (GEOS 627 example)\n",
    "idata = 1\n",
    "\n",
    "X, dlabs, dlabslong, vlabs, vlabslong, ndata, nparm = load_pca_data(idata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9abc688",
   "metadata": {},
   "outputs": [],
   "source": [
    "#==========================================================================\n",
    "### PART 1: MATLAB tutorial\n",
    "# https://www.mathworks.com/help/stats/quality-of-life-in-u-s-cities.html\n",
    "\n",
    "# Load data\n",
    "# idata = 1: City ratings data\n",
    "# idata = 2: Protein consumption data\n",
    "idata = 1\n",
    "X, dlabs, dlabslong, vlabs, vlabslong, ndata, nparm = load_pca_data(idata)\n",
    "\n",
    "# Make a boxplot to look at the distribution of the protein consumption data\n",
    "plt.figure(figsize=(10,8))\n",
    "plt.boxplot(X,vert=False,labels=vlabslong)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35c0b7b9-2cac-40fb-8128-0a66ac36eb4e",
   "metadata": {},
   "outputs": [],
   "source": [
    "C = np.corrcoef(X,rowvar=False)\n",
    "showmat(C,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f5dff05-f7db-4a80-aa7c-33899bbc12d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "s = np.std(X,axis=0,ddof=1)\n",
    "showmat([s],2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c41a4a0-0cd0-4fbe-99e5-bd21157ef774",
   "metadata": {},
   "outputs": [],
   "source": [
    "w = 1/s  # this differs from w in MATLAB code -- we implement weights manually \n",
    "W = np.tile(w,(ndata,1))\n",
    "print(X.shape)\n",
    "print(W.shape)\n",
    "#showmat(W,2)\n",
    "Xprime = X*W    # this is NOT matrix multiplication"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a74db82",
   "metadata": {},
   "outputs": [],
   "source": [
    "pca = PCA()\n",
    "pca.fit(Xprime)\n",
    "coefforth = pca.components_.T  # \n",
    "wcoeff = (np.tile(s.flatten(),(nparm,1))*LA.inv(coefforth)).T  # wcoeff from MATLAB weighted PCA\n",
    "score = pca.fit_transform(Xprime)  # score\n",
    "latent = np.reshape(pca.explained_variance_, (nparm,1))  # latent\n",
    "cscores = zscore(X,ddof=1)@coefforth\n",
    "# how to get tsquared?\n",
    "explained = pca.explained_variance_ratio_ * 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5227813",
   "metadata": {},
   "outputs": [],
   "source": [
    "# visualizing the results\n",
    "plt.figure()\n",
    "plt.plot(score[:,0],score[:,1],'+')\n",
    "plt.xlabel('1st Principal Component')\n",
    "plt.ylabel('2nd Principal Component')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9d3ca37",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.bar(range(1,nparm+1),explained, color='C0')\n",
    "ax.set_xlim([0,10])\n",
    "ax.set_xticks(range(1,nparm+1))\n",
    "ax.tick_params(axis='y', colors='C0')\n",
    "ax.set_xlabel('Principal Component')\n",
    "ax.set_ylabel('Variance Explained (%)',color='C0')\n",
    "ax2 = ax.twinx()\n",
    "ax2.plot(range(0,nparm+1),[0, *list(np.cumsum(explained))], color=\"C3\", marker=\"D\", ms=7)\n",
    "ax2.set_ylim([0,110])\n",
    "ax2.tick_params(axis='y', colors='C3')\n",
    "ax2.set_ylabel('Cumulative Percentage (%)', color='C3')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8bec9e45",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create biplot function similar to MATLAB's\n",
    "def create_biplot(score,coeff,labels=None):\n",
    "    plt.figure(figsize=(10,8))\n",
    "    xs = score[:,0]\n",
    "    ys = score[:,1]\n",
    "    n = coeff.shape[0]\n",
    "    scalex = 1.0/(xs.max() - xs.min())\n",
    "    scaley = 1.0/(ys.max() - ys.min())\n",
    "    plt.axvline(0,color='k',alpha=0.5)\n",
    "    plt.axhline(0,color='k',alpha=0.5)\n",
    "    plt.scatter(xs * scalex,ys * scaley,s=7, color='red')\n",
    "    for i in range(n):\n",
    "        plt.arrow(0, 0, coeff[i,0], coeff[i,1], color='blue', alpha=0.8)\n",
    "        if labels is None:\n",
    "            plt.text(coeff[i,0]*1.05, coeff[i,1]*1.05, \"Var\"+str(i+1), color='k', ha='center', va='center')\n",
    "        else:\n",
    "            plt.text(coeff[i,0]*1.05, coeff[i,1]*1.05, labels[i], color='k', ha='center', va='center')\n",
    "    plt.xlabel(\"Principal Component {}\".format(1))\n",
    "    plt.ylabel(\"Principal Component{}\".format(2))\n",
    "    plt.grid()\n",
    "\n",
    "# Create biplot for first two principal components\n",
    "create_biplot(score[:,0:2],np.transpose(coefforth.T[0:2, :]),vlabslong)\n",
    "plt.xlim([-0.55,0.55])\n",
    "plt.ylim([-0.55,0.55])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f16ebce",
   "metadata": {},
   "outputs": [],
   "source": [
    "#==========================================================================\n",
    "### PART 2\n",
    "\n",
    "# Load data\n",
    "# idata = 1: City ratings data\n",
    "# idata = 2: Protein consumption data\n",
    "idata = 2\n",
    "X, dlabs, dlabslong, vlabs, vlabslong, ndata, nparm = load_pca_data(idata)\n",
    "\n",
    "# Make a boxplot to look at the distribution of the protein consumption data\n",
    "plt.figure(figsize=(10,8))\n",
    "plt.boxplot(X,vert=False,labels=vlabslong)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8b1967b-560a-44fc-bd8c-b251138af3d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "n,p = np.shape(X)\n",
    "xinds = range(1,nparm+1)\n",
    "X0 = X\n",
    "\n",
    "# correlation coefficients\n",
    "R = np.corrcoef(X0,rowvar=False)\n",
    "R2 = np.triu(R)\n",
    "Matrix(R2)\n",
    "\n",
    "# display correlation coefficients '\n",
    "print(''.join(['%10s',*['%8s' for i in range(p)]]) % ('',*vlabs))\n",
    "for ii in range(p):\n",
    "    print(''.join(['%10s',*['%8.3f' for i in range(p)]]) % (vlabslong[ii],*R2[ii,:]))\n",
    "    \n",
    "print('\\nlargest correlation coefficient magnitudes:')\n",
    "Rtemp = R2\n",
    "Rtemp[np.where(Rtemp>0.999)] = np.nan\n",
    "Matrix(Rtemp)\n",
    "for kk in range(20):\n",
    "    val = np.nanmax(abs(Rtemp))\n",
    "    ii,jj = np.unravel_index(np.nanargmax(abs(Rtemp)), Rtemp.shape)\n",
    "    # _,p = stats.pearsonr(X[:,ii],X[:,jj])\n",
    "    _,pval,rlo,rup = pearsonr_ci(X[:,ii],X[:,jj],alpha=0.05)\n",
    "    print('  corr[%10s,%10s] = %6.3f  [%6.3f, %6.3f], p = %6.3f' % \n",
    "          (vlabslong[ii],vlabslong[jj],Rtemp[ii,jj],rlo,rup,pval))\n",
    "    Rtemp[ii,jj] = np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65904429-42fd-4bf2-b4ef-d1217b1101ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "# scatterplot matrix: uncentered, unstandardized\n",
    "df = pd.DataFrame(X0,columns=vlabslong)\n",
    "scatter_matrix(df, alpha=1, figsize=(8,8), diagonal='hist')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "354c02c8-51cb-4812-9835-23d54280cf1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exploring relationships among eig, cov, svd, pca, etc\n",
    "# from wiki: http://en.wikipedia.org/wiki/Principal_component_analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "322ded41-9f2b-4918-9368-3782ddece3a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "### REMOVE THE MEAN FOR EACH COLUMN\n",
    "u = np.mean(X,axis=0)         # dimensionless vector\n",
    "print(vlabslong)\n",
    "showmat([u],2)\n",
    "u = np.reshape(u,(1,len(u)))  # row vector\n",
    "ones = np.ones((ndata,1))        # column vector\n",
    "print(ones.shape)\n",
    "print(u.shape)\n",
    "print((ones@u).shape)\n",
    "#showmat(ones@u,2)\n",
    "# Craft the centered matrix\n",
    "B = X - ones@u\n",
    "#showmat(B,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9dce57ae-b07a-4616-b2cb-844fb868c457",
   "metadata": {},
   "outputs": [],
   "source": [
    "### COVARIANCE MATRIX\n",
    "C = (1/(ndata-1)) * (B.T @ B)  \n",
    "# LA.norm(C-np.cov(B,rowvar=False)) / LA.norm(C) # check\n",
    "# LA.norm(C-np.cov(X,rowvar=False)) / LA.norm(C) # check\n",
    "\n",
    "# standard deviations of each column of X\n",
    "s = np.sqrt(np.diag(C))        # dimensionless vector\n",
    "s = np.reshape(s, (1,len(s)))  # row vector\n",
    "# s = np.sqrt(np.var(X, axis=0, ddof=1))) # alternative\n",
    "showmat(s,2)\n",
    "\n",
    "# important: np.sqrt and np.var differ from MATLAB's sqrt and var \n",
    "# due to default setting of ddof=0\n",
    "\n",
    "# check for sigma1\n",
    "Cdiag = np.diag(np.diag(C))\n",
    "hCdiag = LA.sqrtm(Cdiag)\n",
    "# print([np.std(X[:,0], ddof=1), hCdiag[0][0]]) # check\n",
    "# print([np.var(X[:,0], ddof=1), Cdiag[0][0]]) # check"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0315bea-799d-48bd-bed7-1a62dbb38dea",
   "metadata": {},
   "outputs": [],
   "source": [
    "### CORRELATION MATRIX (is a scaled covariance matrix)\n",
    "# R = np.corrcoef(X, rowvar=False) # alternative\n",
    "# R = np.corrcoef(B, rowvar=False) # alternative\n",
    "# R = corrcov(C) # alternative\n",
    "# R = C/(s.T @ s) # alternative\n",
    "# R = np.diag(1/s.flatten()) @ C @ np.diag(1/s.flatten()) # alternative\n",
    "R = LA.inv(hCdiag) @ C @ LA.inv(hCdiag)\n",
    "# Ccheck = hCdiag @ R @ hCdiag # check\n",
    "# LA.norm(Ccheck - C) / LA.norm(C) # check"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "077867d5-950b-4ce7-b904-f65b80ec6196",
   "metadata": {},
   "outputs": [],
   "source": [
    "### STANDARDIZED (AND CENTERED) MATRIX\n",
    "# Z = B / (h@s) # alternative\n",
    "# Z = B @ np.diag(1/s.flatten()) # alternative\n",
    "Z = B @ LA.inv(hCdiag)\n",
    "# LA.norm(Z - zscore(X,ddof=1)) / LA.norm(Z) # check\n",
    "# Bcheck = Z * (h@s) # check reverse\n",
    "# LA.norm(B-Bcheck) / LA.norm(B) # check reverse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8dcd1980",
   "metadata": {},
   "outputs": [],
   "source": [
    "### EIGENVALUE DECOMPOSITION\n",
    "eigvalC, Vc = LA.eig(C)\n",
    "Dc = np.real(np.diag(eigvalC))\n",
    "# LA.norm(C@Vc - Vc@Dc) / LA.norm(C@Vc) # check\n",
    "## if needed: sort eigenvalues and rearrange V\n",
    "isort = np.argsort(eigvalC)[::-1] # [::-1] flips the array\n",
    "Vc = Vc[:,isort]\n",
    "Dc = np.diag(eigvalC[isort])\n",
    "# LA.norm(C@Vc - Vc@Dc) / LA.norm(C@Vc) # check\n",
    "eigvalC = eigvalC[isort]  # dimensionless\n",
    "eigvalC = np.reshape(eigvalC,(nparm,1))  # column vector\n",
    "showmat(eigvalC.T,2)\n",
    "\n",
    "eigvalR, Vr = LA.eig(R)\n",
    "Dr = np.real(np.diag(eigvalR))\n",
    "isort = np.argsort(eigvalR)[::-1] # [::-1] flips the array\n",
    "Vr = Vr[:,isort]\n",
    "Dr = np.diag(eigvalR[isort])\n",
    "# LA.norm(R@Vr - Vr@Dr) / LA.norm(R@Vr) # check\n",
    "eigvalR = eigvalR[isort]  # dimensionless\n",
    "eigvalR = np.reshape(eigvalR,(nparm,1))  # column vector\n",
    "showmat(eigvalR.T,2)\n",
    "# Vr @ Dr @ Vr.T\n",
    "# LA.inv(hCdiag) @ Vc @ Dc @ (LA.inv(hCdiag) @ Vc).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca5c431b",
   "metadata": {},
   "outputs": [],
   "source": [
    "experiment = False\n",
    "\n",
    "if experiment:\n",
    "    # experiment with multiplying each column of A by a factor\n",
    "    A = np.random.rand(8,3)\n",
    "    s = np.sqrt(np.var(A, axis=0, ddof=1))  # dimensionless\n",
    "    s = np.reshape(s,(1,3))  # row vector\n",
    "    h = np.ones((8,1))  # column vector\n",
    "    print('A = \\n', A)\n",
    "    # this looks like dividing each column by its standard deviation\n",
    "    ans = A / (h@s)\n",
    "    print('A / (h@s) = \\n', ans)\n",
    "    ans = A / np.tile(s,(8,1))\n",
    "    print('A / np.tile(s,(8,1)) = \\n', ans)\n",
    "    # this is the series of elemetary matrix operations\n",
    "    ans = A @ np.diag([1/s[0,0],1,1]) @ np.diag([1,1/s[0,1],1]) @ np.diag([1,1,1/s[0,2]])\n",
    "    # this is what it means\n",
    "    ans = A @ np.diag(1/s.flatten())\n",
    "    print('A @ np.diag(1/s.flatten()) = \\n', ans)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0d6049f",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('mean(X) =')\n",
    "showmat([np.mean(X,axis=0)],2)\n",
    "print('mean(B) =')\n",
    "showmat([np.mean(B,axis=0)],2)\n",
    "print('mean(Z) =')\n",
    "showmat([np.mean(Z,axis=0)],2)\n",
    "print('var(X) =')\n",
    "showmat([np.var(X,axis=0,ddof=1)],2)\n",
    "print('var(B) =')\n",
    "showmat([np.var(B,axis=0,ddof=1)],2)\n",
    "print('var(Z) =')\n",
    "showmat([np.var(Z,axis=0,ddof=1)],2)\n",
    "print('std(X) =')\n",
    "showmat([np.std(X,axis=0,ddof=1)],2)\n",
    "print('std(B) =')\n",
    "showmat([np.std(B,axis=0,ddof=1)],2)\n",
    "print('std(Z) =')\n",
    "showmat([np.std(Z,axis=0,ddof=1)],2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8ede487",
   "metadata": {},
   "outputs": [],
   "source": [
    "### SINGULAR VALUE DECOMPOSITION of B\n",
    "# the scores matrix is nothing more than U@S from the SVD\n",
    "Ub,Sb,Vb = svdmat(B)\n",
    "svalb = np.diag(Sb)  # dimensionless\n",
    "svalb = np.reshape(svalb,(nparm,1))  # column vector\n",
    "USb = Ub@Sb\n",
    "# check singular values with eigenvalues of covariance matrix\n",
    "print(LA.norm(np.square(svalb)/(ndata-1) - eigvalC) / LA.norm(eigvalC))\n",
    "# check: Vb = Vc. Use abs() to allow for sign flips\n",
    "# abs(Vc) - abs(Vb)\n",
    "\n",
    "### SINGULAR VALUE DECOMPOSITION of Z\n",
    "# Vz = Vr (allowing sign changes)\n",
    "Uz,Sz,Vz = svdmat(Z)\n",
    "svalz = np.diag(Sz)  # dimensionless\n",
    "svalz = np.reshape(svalz,(nparm,1))  # column vector\n",
    "USz = Uz@Sz\n",
    "# check singular values with eigenvalues of correlation matrix\n",
    "print(LA.norm(np.square(svalz)/(ndata-1) - eigvalR) / LA.norm(eigvalR))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1952bbf9",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Test 1: Use centered matrix\n",
    "# VB = Vb = Vc (allowing for some sign flips on columns of V)\n",
    "# USB = Ub@Sb\n",
    "# pcvarB = eigvalC\n",
    "pca = PCA()\n",
    "pca.fit(B)\n",
    "VB = pca.components_.T  # coeff\n",
    "USB = pca.fit_transform(B)  # score\n",
    "pcvarB = np.reshape(pca.explained_variance_, (nparm,1))  # latent\n",
    "showmat(pcvarB.T,2)\n",
    "showmat(eigvalC.T,2)\n",
    "# this is equivalent, since pca will center the matrix (i.e., remove mean)\n",
    "# pca = PCA()\n",
    "# pca.fit(X)\n",
    "# V = pca.components_.T  # coeff\n",
    "# US = pca.fit_transform(X)  # score\n",
    "\n",
    "# Bcheck = USB @ VB.T\n",
    "# LA.norm(B-Bcheck) / LA.norm(B)\n",
    "# # orthonormal\n",
    "# LA.norm(VB.T @ VB - np.eye(nparm))\n",
    "# US1_check = B @ VB\n",
    "# LA.norm(USB - US1_check) / LA.norm(USB)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57229d59",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Test 2 (example used in matlab)\n",
    "# note: this gives different US and V from Test 1\n",
    "w = 1/s  # this differs from w in MATLAB code -- we implement weights manually\n",
    "W = np.tile(w,(ndata,1))\n",
    "Bprime = B*W\n",
    "pca = PCA()\n",
    "pca.fit(Bprime)\n",
    "Vw_not_orthonormal = pca.components_.T  # \n",
    "Vw = (np.tile(s.flatten(),(nparm,1))*LA.inv(Vw_not_orthonormal)).T  # wcoeff from MATLAB weighted PCA\n",
    "USw = pca.fit_transform(Bprime)  # score\n",
    "pcvarBw = np.reshape(pca.explained_variance_, (nparm,1))  # latent\n",
    "\n",
    "# % this is equivalent\n",
    "# %[Vw,USw] = pca(B,'VariableWeights','variance')\n",
    "# Bcheck = USw * Vw';\n",
    "# norm(B - Bcheck) / norm(B)\n",
    "# disp('Vw is NOT orthornomal:');\n",
    "# norm(Vw'*Vw - eye(nparm))\n",
    "# Vworth = inv(hCdiag)*Vw;   % note: inv(hCdiag) = diag(1./s)\n",
    "# % orthonormal:\n",
    "# norm(Vworth'*Vworth - eye(nparm))\n",
    "# % this shows how USw can be computed\n",
    "# USw_check = Z*Vworth;\n",
    "# norm(USw - USw_check) / norm(USw)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34461b39",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Test 3: Use centered+standardized matrix as input\n",
    "# VZ = Vz = Vr = LA.inv(hCdiag)@Vw  (allowing for sign flips)\n",
    "# USZ = USw (allowing for sign flips)\n",
    "# pcvarZ = eigvalR\n",
    "pca = PCA()\n",
    "pca.fit(Z)\n",
    "VZ = pca.components_.T  # coeff\n",
    "USZ = pca.fit_transform(Z)  # score\n",
    "pcvarZ = np.reshape(pca.explained_variance_, (nparm,1))  # latent\n",
    "\n",
    "# Zcheck = USZ @ VZ.T\n",
    "# LA.norm(Z-Zcheck) / LA.norm(Z)\n",
    "# Bcheck = USZ @ VZ.T @ hCdiag\n",
    "# LA.norm(B-Bcheck) / LA.norm(B)\n",
    "# # orthonormal\n",
    "# LA.norm(VZ.T @ VZ - np.eye(nparm))\n",
    "# VZ_check = LA.inv(hCdiag)@Vw\n",
    "# LA.norm(VZ - VZ_check) / VZ.norm(USB)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "71dbba44",
   "metadata": {},
   "outputs": [],
   "source": [
    "#==============================\n",
    "# tips for the lab exercise (lab_pca.m)\n",
    "\n",
    "# cumulative variance\n",
    "for kk in range(2):\n",
    "    if kk == 0:\n",
    "        pcvar = pcvarB\n",
    "        stlab = 'centered'\n",
    "    else:\n",
    "        pcvar = pcvarZ\n",
    "        stlab = 'centered+standardized'\n",
    "    propvar = pcvar/np.sum(pcvar)\n",
    "    # from the matlab tutorial: explained = 100*cpropvar (idata=1)\n",
    "    cpropvar = np.cumsum(propvar)\n",
    "    print('\\n')\n",
    "    print('Importance of principal components [%s]' % stlab)\n",
    "    print('Std-Dev  : sqrt( lam_k )')\n",
    "    print('Var      : lam_k = s_k^2 / (ndata-1) : eigenvalues of the covariance (or correlation) matrix of X')\n",
    "    print('Prop-Var : proportion of variance')\n",
    "    print('Std-Dev  : cumulative proportion of variance')\n",
    "    print('----------------------------------------------')\n",
    "    print('PC#\\tStd-Dev\\t\\tVar\\t\\tProp-Var\\tCum-Prop')\n",
    "    for i in range(nparm):\n",
    "        print('%d\\t%f\\t%f\\t%f\\t%f' % (i+1,np.sqrt(pcvar[i]),pcvar[i],propvar[i],cpropvar[i]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8a834d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot components\n",
    "fig=plt.figure(figsize=(10,14))\n",
    "for ii in range(nparm):\n",
    "    plt.subplot(5,2,ii+1)\n",
    "    plt.plot([0.5,nparm+0.5],[0,0],'k')\n",
    "    h1 = plt.plot(range(1,nparm+1),VB[:,ii],'r.-',label='centered')\n",
    "    h2 = plt.plot(range(1,nparm+1),VZ[:,ii],'b.-',label='standardized')\n",
    "    plt.axis([0,nparm+1,-1,1])\n",
    "    plt.title('PC-%d' % (ii+1))\n",
    "    plt.xticks(range(1,nparm+1),labels=vlabs)\n",
    "    plt.grid()\n",
    "    if ii == 0:\n",
    "        plt.legend()\n",
    "        plt.subplots_adjust(wspace=0.2,hspace=0.5)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "inverse [conda env:inverse]",
   "language": "python",
   "name": "conda-env-inverse-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}