.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "api/gallery/statistics/pca/plot_01_pearson_hotelling_pca.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code or to run this example in your browser via JupyterLite. .. rst-class:: sphx-glr-example-title .. _sphx_glr_api_gallery_statistics_pca_plot_01_pearson_hotelling_pca.py: Pearson and Hotelling's principal component analysis ===================================================== Pearson (1901) asked for the line that best fits a cloud of points when both coordinates carry error: the one minimizing the *perpendicular* distances, not the vertical ones of least-squares regression. It runs along the leading eigenvector of the covariance matrix. Hotelling (1933) turned the same eigenvectors into a method for many variables: replace correlated test scores by a few uncorrelated "components" that carry most of their variance. .. GENERATED FROM PYTHON SOURCE LINES 15-22 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np from mathematicskit.statistics import linear_regression, principal_component_analysis rng = np.random.default_rng(0) .. GENERATED FROM PYTHON SOURCE LINES 23-29 Pearson's line of closest fit ------------------------------- Regressing :math:`y` on :math:`x` and :math:`x` on :math:`y` give two different lines; the first principal axis lies between them and minimizes the sum of squared perpendicular distances. .. GENERATED FROM PYTHON SOURCE LINES 29-62 .. code-block:: Python latent = rng.normal(size=300) points = np.column_stack([latent + 0.6 * rng.normal(size=300), 0.8 * latent + 0.6 * rng.normal(size=300)]) pca = principal_component_analysis(points) direction, center = pca.components[0], pca.mean y_on_x = linear_regression(points[:, 0], points[:, 1]).coefficients x_on_y = linear_regression(points[:, 1], points[:, 0]).coefficients def perpendicular_ss(slope, intercept): return np.sum((points[:, 1] - slope * points[:, 0] - intercept) ** 2) / (1 + slope**2) pca_slope = direction[1] / direction[0] lines = { "y on x (OLS)": (y_on_x[1], y_on_x[0]), "x on y (OLS)": (1 / x_on_y[1], -x_on_y[0] / x_on_y[1]), "Pearson / first principal axis": (pca_slope, center[1] - pca_slope * center[0]), } for name, (slope, intercept) in lines.items(): print(f"{name:32s} slope {slope:5.3f}, perpendicular sum of squares {perpendicular_ss(slope, intercept):7.2f}") print(f"smallest eigenvalue x (n - 1) = {pca.explained_variance[1] * (len(points) - 1):7.2f}") fig, (ax0, ax1, ax2) = plt.subplots(1, 3, figsize=(15, 4.5)) ax0.scatter(*points.T, s=6, color="0.6") xs = np.linspace(points[:, 0].min(), points[:, 0].max(), 2) for (name, (slope, intercept)), style in zip(lines.items(), ["--", ":", "-"], strict=True): ax0.plot(xs, slope * xs + intercept, style, lw=2, label=name) ax0.set_aspect("equal") ax0.set_title("Pearson (1901): the line of closest fit") ax0.legend(fontsize=8) .. image-sg:: /api/gallery/statistics/pca/images/sphx_glr_plot_01_pearson_hotelling_pca_001.png :alt: Pearson (1901): the line of closest fit :srcset: /api/gallery/statistics/pca/images/sphx_glr_plot_01_pearson_hotelling_pca_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none y on x (OLS) slope 0.627, perpendicular sum of squares 97.96 x on y (OLS) slope 1.135, perpendicular sum of squares 108.00 Pearson / first principal axis slope 0.796, perpendicular sum of squares 91.04 smallest eigenvalue x (n - 1) = 91.04 .. GENERATED FROM PYTHON SOURCE LINES 63-70 Hotelling's components of six test scores ------------------------------------------- Six scores driven by two hidden abilities, "verbal" and "quantitative", plus noise. Standardizing (the correlation matrix) puts the tests on an equal footing. Two components carry most of the variance, and their loadings separate the two kinds of test. .. GENERATED FROM PYTHON SOURCE LINES 70-99 .. code-block:: Python n = 500 verbal, quantitative = rng.normal(size=(2, n)) loadings = np.array([[0.9, 0.1], [0.8, 0.2], [0.85, 0.0], [0.1, 0.9], [0.2, 0.8], [0.0, 0.85]]) scores = np.column_stack([verbal, quantitative]) @ loadings.T + 0.4 * rng.normal(size=(n, 6)) tests = ["reading", "vocabulary", "grammar", "arithmetic", "algebra", "geometry"] hotelling = principal_component_analysis(scores, standardize=True) print("explained variance ratio:", hotelling.explained_variance_ratio.round(3)) ax1.bar(np.arange(1, 7), hotelling.explained_variance_ratio, color="tab:blue") ax1.plot(np.arange(1, 7), np.cumsum(hotelling.explained_variance_ratio), "ko-", label="cumulative") ax1.set_xlabel("component") ax1.set_ylabel("fraction of variance") ax1.set_title("Scree plot: two components dominate") ax1.legend() first, second = hotelling.components[0], hotelling.components[1] ax2.axhline(0, color="0.8") ax2.axvline(0, color="0.8") for name, a, c in zip(tests, first, second, strict=True): ax2.arrow(0, 0, a, c, head_width=0.02, color="tab:red" if name in tests[:3] else "tab:green", length_includes_head=True) ax2.annotate(name, (a, c), textcoords="offset points", xytext=(4, 2), fontsize=8) ax2.set_xlabel("loading on component 1 (general ability)") ax2.set_ylabel("loading on component 2 (verbal vs. quantitative)") ax2.set_title("Hotelling (1933): component loadings") ax2.set_aspect("equal") fig.tight_layout() plt.show() .. rst-class:: sphx-glr-script-out .. code-block:: none explained variance ratio: [0.539 0.348 0.032 0.031 0.028 0.023] .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.084 seconds) .. _sphx_glr_download_api_gallery_statistics_pca_plot_01_pearson_hotelling_pca.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: lite-badge .. image:: images/jupyterlite_badge_logo.svg :target: ../../../../lite/lab/index.html?path=api/gallery/statistics/pca/plot_01_pearson_hotelling_pca.ipynb :alt: Launch JupyterLite :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_01_pearson_hotelling_pca.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_01_pearson_hotelling_pca.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_01_pearson_hotelling_pca.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_