icaoberg / Statistics - Principal Component Analysis (PCA)

Created Sun, 17 May 2026 00:00:00 +0000 Modified Sun, 17 May 2026 17:14:33 -0400

Jupyter Notebook python

statistics pca dimensionality-reduction scikit-learn iris

Statistics - Principal Component Analysis (PCA)¶

Principal Component Analysis (PCA) is one of the most widely used techniques in data science for dimensionality reduction, visualization, and exploratory data analysis.

This notebook covers:

Definition and background
When and why PCA is useful
A real-life example using the Iris dataset
Step-by-step PCA from scratch and with scikit-learn
Interpreting results

Definition¶

According to the NIST/SEMATECH e-Handbook of Statistical Methods:

"Principal components analysis (PCA) is a multivariate technique that analyzes a data table in which observations are described by several inter-correlated quantitative dependent variables. Its goal is to extract the important information from the table, to represent it as a set of new orthogonal variables called principal components, and to display the pattern of similarity of the observations and of the variables as points in maps."

Source: NIST/SEMATECH e-Handbook of Statistical Methods, Section 6.5.5.2 — Principal Components

In simpler terms, PCA finds the directions (principal components) in which the data varies the most. By projecting data onto these directions, we can represent high-dimensional data in fewer dimensions while retaining as much variance (information) as possible.

Key properties of principal components:¶

They are orthogonal (uncorrelated) to each other.
The first principal component captures the most variance in the data.
Each subsequent component captures the most remaining variance, orthogonal to the previous ones.
The number of components is at most equal to the number of original features.

When Is PCA Useful?¶

PCA is particularly valuable under the following conditions:

Condition	Why PCA Helps
Many features (high dimensionality)	Reduces features while preserving variance — combats the curse of dimensionality
Correlated features	Removes redundancy by transforming to uncorrelated components
Visualization	Projects data to 2D or 3D for human-readable scatter plots
Noise reduction	Low-variance components often capture noise; dropping them denoises the data
Pre-processing for ML models	Speeds up training and can improve generalization when features are correlated

PCA is not ideal when:

Features are already uncorrelated.
Interpretability of individual features must be preserved.
The relationship between features and the target is non-linear (consider kernel PCA or autoencoders instead).

Setup¶

We use standard scientific Python libraries available in most conda/pip environments.

In [1]:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.datasets import load_iris
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA

# Consistent styling
sns.set_theme(style='whitegrid')
plt.rcParams['figure.dpi'] = 120

print('Libraries loaded successfully.')

Libraries loaded successfully.

Real-Life Example: The Iris Dataset¶

The Iris dataset (Fisher, 1936) is a classic benchmark in statistics and machine learning. It describes 150 iris flowers across three species (Iris setosa, Iris versicolor, Iris virginica), measured along four features:

Sepal length (cm)
Sepal width (cm)
Petal length (cm)
Petal width (cm)

Why is this a good PCA example? The four measurements are correlated (petal length and petal width in particular are highly correlated), so PCA can compress the 4-dimensional space into 2 or 3 principal components without losing much information — making the three species visually separable in a 2D plot.

In [2]:

iris = load_iris(as_frame=True)
df = iris.frame
df.columns = ['sepal_length', 'sepal_width', 'petal_length', 'petal_width', 'species']
df['species_name'] = df['species'].map(dict(enumerate(iris.target_names)))

print(f'Shape: {df.shape}')
df.head()

Shape: (150, 6)

Out[2]:

	sepal_length	sepal_width	petal_length	petal_width	species_name
0	5.1	3.5	1.4	0.2	setosa
1	4.9	3.0	1.4	0.2	setosa
2	4.7	3.2	1.3	0.2	setosa
3	4.6	3.1	1.5	0.2	setosa
4	5.0	3.6	1.4	0.2	setosa

Step 1 — Check Feature Correlations¶

PCA is most useful when features are correlated. Let's verify that is the case here.

In [3]:

fig, ax = plt.subplots(figsize=(6, 5))
corr = df[['sepal_length', 'sepal_width', 'petal_length', 'petal_width']].corr()
sns.heatmap(corr, annot=True, fmt='.2f', cmap='coolwarm', center=0, ax=ax)
ax.set_title('Feature Correlation Matrix — Iris Dataset')
plt.tight_layout()
plt.show()

print('\nObservation: petal_length and petal_width are highly correlated (r ≈ 0.96).')
print('sepal_length is also moderately correlated with petal measurements.')
print('This redundancy makes PCA valuable — we can reduce 4 features to 2 components.')

No description has been provided for this image

Observation: petal_length and petal_width are highly correlated (r ≈ 0.96).
sepal_length is also moderately correlated with petal measurements.
This redundancy makes PCA valuable — we can reduce 4 features to 2 components.

Step 2 — Standardize the Features¶

PCA is sensitive to the scale of features. Standardizing (zero mean, unit variance) ensures each feature contributes equally before computing principal components.

In [4]:

features = ['sepal_length', 'sepal_width', 'petal_length', 'petal_width']
X = df[features].values

scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

print('Before scaling:')
print(f'  Mean: {X.mean(axis=0).round(2)}')
print(f'  Std:  {X.std(axis=0).round(2)}')
print('\nAfter scaling:')
print(f'  Mean: {X_scaled.mean(axis=0).round(10)}  (≈ 0)')
print(f'  Std:  {X_scaled.std(axis=0).round(2)}')

Before scaling:
  Mean: [5.84 3.06 3.76 1.2 ]
  Std:  [0.83 0.43 1.76 0.76]

After scaling:
  Mean: [-0. -0. -0. -0.]  (≈ 0)
  Std:  [1. 1. 1. 1.]

Step 3 — Apply PCA and Examine Explained Variance¶

We first fit PCA retaining all components to see how much variance each one explains.

In [5]:

pca_full = PCA()
pca_full.fit(X_scaled)

explained = pca_full.explained_variance_ratio_
cumulative = np.cumsum(explained)

print('Explained variance per component:')
for i, (ev, cv) in enumerate(zip(explained, cumulative), 1):
    print(f'  PC{i}: {ev:.1%}  (cumulative: {cv:.1%})')

# Scree plot
fig, axes = plt.subplots(1, 2, figsize=(11, 4))

axes[0].bar(range(1, len(explained)+1), explained * 100, color='steelblue')
axes[0].set_xlabel('Principal Component')
axes[0].set_ylabel('Explained Variance (%)')
axes[0].set_title('Scree Plot')
axes[0].set_xticks(range(1, len(explained)+1))

axes[1].plot(range(1, len(cumulative)+1), cumulative * 100, marker='o', color='darkorange')
axes[1].axhline(95, linestyle='--', color='gray', label='95% threshold')
axes[1].set_xlabel('Number of Components')
axes[1].set_ylabel('Cumulative Explained Variance (%)')
axes[1].set_title('Cumulative Explained Variance')
axes[1].set_xticks(range(1, len(cumulative)+1))
axes[1].legend()

plt.tight_layout()
plt.show()

Explained variance per component:
  PC1: 73.0%  (cumulative: 73.0%)
  PC2: 22.9%  (cumulative: 95.8%)
  PC3: 3.7%  (cumulative: 99.5%)
  PC4: 0.5%  (cumulative: 100.0%)

Step 4 — Project to 2D and Visualize¶

The first two principal components together explain the majority of the variance. Plotting samples in this 2D space often reveals natural clusters.

In [6]:

pca2 = PCA(n_components=2)
X_pca = pca2.fit_transform(X_scaled)

pca_df = pd.DataFrame(X_pca, columns=['PC1', 'PC2'])
pca_df['species'] = df['species_name'].values

fig, ax = plt.subplots(figsize=(8, 6))
palette = {'setosa': '#e41a1c', 'versicolor': '#377eb8', 'virginica': '#4daf4a'}

for species, group in pca_df.groupby('species'):
    ax.scatter(group['PC1'], group['PC2'], label=species,
               color=palette[species], alpha=0.8, edgecolors='white', s=80)

ev = pca2.explained_variance_ratio_
ax.set_xlabel(f'PC1 ({ev[0]:.1%} variance explained)')
ax.set_ylabel(f'PC2 ({ev[1]:.1%} variance explained)')
ax.set_title('Iris Dataset — 2D PCA Projection')
ax.legend(title='Species')
plt.tight_layout()
plt.show()

print(f'Total variance explained by 2 components: {ev.sum():.1%}')

Total variance explained by 2 components: 95.8%

Step 5 — Interpret the Principal Components (Loadings)¶

The loadings (eigenvectors) tell us how much each original feature contributes to each principal component. Large absolute values indicate strong contribution.

In [7]:

loadings = pd.DataFrame(
    pca2.components_.T,
    index=features,
    columns=['PC1', 'PC2']
)
print('Component Loadings:')
print(loadings.round(3))

fig, ax = plt.subplots(figsize=(7, 4))
loadings.plot(kind='bar', ax=ax, color=['steelblue', 'darkorange'])
ax.axhline(0, color='black', linewidth=0.8)
ax.set_title('PCA Loadings — Contribution of Each Feature')
ax.set_xlabel('Original Feature')
ax.set_ylabel('Loading')
ax.set_xticklabels(features, rotation=30, ha='right')
ax.legend()
plt.tight_layout()
plt.show()

print('\nInterpretation:')
print('  PC1 is driven mainly by petal features (+ sepal_length).')
print('  PC2 captures mostly sepal_width variation.')

Component Loadings:
                PC1    PC2
sepal_length  0.521  0.377
sepal_width  -0.269  0.923
petal_length  0.580  0.024
petal_width   0.565  0.067

Interpretation:
  PC1 is driven mainly by petal features (+ sepal_length).
  PC2 captures mostly sepal_width variation.

Summary¶

Step	Action	Result
1	Checked correlations	Found strong correlation between petal features
2	Standardized data	Equal contribution from all features
3	Fit full PCA	PC1 + PC2 explain ~95% of variance
4	2D projection	Three species visually separable
5	Inspected loadings	PC1 ≈ petal size; PC2 ≈ sepal width

PCA was effective here because:

The original four features are highly correlated.
Two components capture nearly all variance — a 4→2 compression with minimal loss.
The projection reveals meaningful structure (species clusters) invisible in any single feature alone.