Data Analytics with Machine Learning

Scales of Measurement

Before any machine learning model can be constructed, the data must be understood in terms of its scales of measurement. These scales determine which mathematical operations are meaningful and guide feature encoding, distance metrics, and ultimately the choice of algorithms. There are four fundamental scales.

Nominal attributes are categorical labels without any intrinsic ordering. Examples include country of origin, product colour, or a unique identifier. The only permissible operation is equality testing ($x_i = x_j$). For machine learning, nominal variables are typically transformed using one-hot encoding, which creates $k$ binary dummy variables for $k$ distinct categories, ensuring that no artificial ordinal relationship is imposed. Note that one-hot encoding increases dimensionality and, if categories are numerous, may require dimensionality reduction or target encoding.

Ordinal attributes possess a meaningful order but the intervals between successive values are not necessarily uniform. A credit rating from AAA to C, or a Likert scale from "strongly disagree" to "strongly agree" are ordinal. While we can state that AAA $\succ$ AA $\succ$ A, we cannot claim that the difference between AAA and AA equals the difference between B and CCC. Ordinal encoding, where integers preserve the order (e.g., 1, 2, 3, …), is common, but caution is required when using distance-based models, which implicitly treat the encoded numbers as interval-scaled.

Interval-scaled attributes have equal and meaningful differences, but no true zero point. Temperature in Celsius or Fahrenheit is interval: the difference between $20^\circ$C and $30^\circ$C is the same as between $10^\circ$C and $20^\circ$C, but $0^\circ$C does not imply the absence of heat. All arithmetic operations (addition, subtraction) are valid, but ratios are not. Standardization (z-score $X' = \frac{X - \mu}{\sigma}$) is applied to interval variables to remove scale effects.

Ratio-scaled attributes have all the properties of interval data plus a true zero, allowing meaningful ratios. Examples include income, height, elapsed time, and quantity sold. A household earning \100,000$earns twice as much as one earning$$50,000$. Ratio variables are also amenable to standardization, but logarithmic or Box-Cox transformations are often applied to handle positive skewness typical in economic indicators.

The hierarchy of measurement scales — nominal $\subset$ ordinal $\subset$ interval $\subset$ ratio — dictates permissible transformations. Nominal data allows only one-to-one mappings; ordinal data permits monotonic increasing transformations; interval data supports linear transformations $y = ax + b$ where $a > 0$; and ratio data allows similarity transformations $y = ax$ where $a > 0$.

Feature Engineering

Feature engineering is the process of transforming raw data into a representation that better captures the underlying structure of the problem, thereby improving model accuracy, interpretability, and generalizability. Formally, given a dataset $\mathcal{D} = \{(\mathbf{x}_i, y_i)\}_{i=1}^n$ with original feature vectors $\mathbf{x}_i \in \mathbb{R}^d$, feature engineering constructs a mapping $\phi: \mathbb{R}^d \to \mathbb{R}^m$ such that a learning algorithm $f: \mathbb{R}^m \to \mathbb{R}$ achieves lower empirical risk.

Essential Steps in Feature Engineering

1. Exploratory Data Analysis (EDA) precedes any transformation. EDA employs summary statistics, correlation matrices, and visualizations (histograms, scatter plots, box plots) to reveal distributions, outliers, and relationships.

   Univariate analysis examines individual variables through measures of central tendency ($\mu = \frac{1}{n}\sum x_i$, median, mode) and dispersion ($\sigma^2 = \frac{1}{n}\sum (x_i - \mu)^2$, standard deviation, IQR). Histograms, box plots, and kernel density estimates reveal distributional properties such as skewness $\gamma_1 = \frac{E[(X-\mu)^3]}{\sigma^3}$ and kurtosis $\gamma_2 = \frac{E[(X-\mu)^4]}{\sigma^4} - 3$.

   Bivariate analysis explores relationships between variables. The Pearson correlation coefficient $r_{xy} = \frac{\sum (x_i - \bar{x})(y_i - \bar{y})}{\sqrt{\sum (x_i - \bar{x})^2}\sqrt{\sum (y_i - \bar{y})^2}}$ quantifies linear associations. However, correlation does not imply causation, and nonlinear relationships may exhibit $r \approx 0$ despite strong dependencies.

   Missing data analysis is critical: the pattern of missingness — whether missing completely at random (MCAR), missing at random (MAR), or missing not at random (MNAR) — determines appropriate imputation strategies. Simple mean imputation $\hat{x}_i = \bar{x}$ preserves the sample mean but artificially reduces variance. More sophisticated approaches include k-nearest neighbors imputation, multiple imputation by chained equations (MICE), or model-based methods.

   Outlier detection employs both statistical and visual techniques. The IQR method flags observations outside $[Q_1 - 1.5 \times IQR, Q_3 + 1.5 \times IQR]$ as potential outliers. Z-score methods identify points where $|z_i| = \left|\frac{x_i - \mu}{\sigma}\right| > 3$.

2. Handling Categorical Variables is perhaps the most common transformation. Nominal variables are encoded via one-hot encoding, which for a categorical feature with $k$ levels yields a binary indicator matrix $\mathbf{I} \in \{0,1\}^{n \times k}$. For high-cardinality features, label encoding (assigning integers $0,1,\dots,k-1$) may be used with tree-based models. An intermediate approach is target encoding, where each category is replaced by the mean of the target variable for that category, often with smoothing to avoid overfitting.

3. Numerical Transformations and Scaling are essential for distance-based algorithms like $k$-NN, SVM with RBF kernel, and neural networks:

   • Standardization (z-score): $x'_i = \frac{x_i - \mu}{\sigma}$, yielding zero mean and unit variance.
   • Min-max scaling: $x'_i = \frac{x_i - \min(X)}{\max(X) - \min(X)}$, bounding values to $[0,1]$.
   • Robust scaling: $x'_i = \frac{x_i - \text{median}}{\text{IQR}}$, providing resilience against extreme values.
   • Logarithmic or Box-Cox transformations to reduce skewness and stabilize variance: $x'(\lambda) =$ \begin{cases}\frac{x^\lambda - 1}{\lambda} & \text{if } \lambda \neq 0 \ \log(x) & \text{if } \lambda = 0\end{cases} $$.

4. Feature Creation and Selection go beyond encoding. Polynomial features $[x_1, x_2, x_1^2, x_1 x_2, x_2^2]$ capture interactions and non-linearities. After creation, feature selection techniques such as recursive feature elimination, regularization ($L_1$ penalty), and mutual information $I(X;Y) = \sum_{y \in Y}\sum_{x \in X}p(x,y)\log\left(\frac{p(x,y)}{p(x)p(y)}\right)$ prune irrelevant or redundant variables.

Performance Evaluation Metrics

Regression Metrics

Given true values $y_i$ and predicted values $\hat{y}_i$, the most fundamental error measure is the mean squared error (MSE):

$$\text{MSE} = \frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2.$$

Because MSE squares errors, it heavily penalizes large deviations. Taking the square root yields the root mean squared error (RMSE):

$$\text{RMSE} = \sqrt{\text{MSE}} = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2}.$$

RMSE is expressed in the same units as the target variable, facilitating interpretation. The quadratic nature makes RMSE differentiable and convex, which is why it is the default optimization objective for many algorithms. However, the sensitivity to outliers can be problematic when predicting variables like income, where a few extreme values may dominate the error. In such cases, the mean absolute error (MAE),

$$\text{MAE} = \frac{1}{n} \sum_{i=1}^{n} |y_i - \hat{y}_i|,$$

provides a more robust alternative because it weights all errors linearly. The coefficient of determination ($R^2$),

$$R^2 = 1 - \frac{\sum_{i=1}^{n} (y_i - \hat{y}_i)^2}{\sum_{i=1}^{n} (y_i - \bar{y})^2},$$

measures the proportion of variance explained, with $R^2=1$ indicating perfect fit. It is scale-free, allowing comparisons across problems.

Classification Metrics

For binary classification, the confusion matrix is the foundation. Each prediction falls into one of four categories: true positive (TP), false positive (FP), false negative (FN), and true negative (TN). From these we define:

• Accuracy: $\frac{TP + TN}{TP + FP + FN + TN}$. Misleading when classes are imbalanced.
• Precision (positive predictive value): $\frac{TP}{TP + FP}$. Measures how many of the predicted positives are actually positive. High precision is crucial when the cost of a false positive is high.
• Recall (sensitivity, true positive rate): $\frac{TP}{TP + FN}$. Captures how many of the actual positives are correctly identified. In medical diagnosis, high recall minimizes the risk of missing a disease.
• Specificity (true negative rate): $\frac{TN}{TN + FP}$. The analogue of recall for the negative class.
• F1-Score: $2 \cdot \frac{\text{Precision} \cdot \text{Recall}}{\text{Precision} + \text{Recall}}$. The harmonic mean of precision and recall.

In economic terms, the choice between precision and recall maps to the relative costs of Type I errors (false positives) and Type II errors (false negatives). A spam filter with high precision avoids blocking legitimate emails (low false positives); a cancer screening with high recall avoids missing cases (low false negatives).

Advanced Classification Metrics

The Receiver Operating Characteristic (ROC) curve plots the true positive rate (recall) against the false positive rate ($1 - \text{specificity}$) at various threshold settings. The Area Under the ROC Curve (AUC-ROC) is a threshold-invariant metric; an AUC of 1.0 indicates perfect discrimination, while 0.5 indicates random guessing. The Precision-Recall (PR) curve is preferred when classes are highly imbalanced.

Log Loss (Cross-Entropy) for probabilistic classifiers:

$$\text{LogLoss} = -\frac{1}{n}\sum_{i=1}^n \left[y_i \log(\hat{p}_i) + (1-y_i)\log(1-\hat{p}_i)\right],$$

where $\hat{p}_i$ is the predicted probability. Log loss heavily penalizes confident but wrong predictions, making it appropriate when probabilistic calibration matters.

Holdout and Cross-Validation

Holdout validation splits data into training (typically 70-80%) and test sets. The k-fold cross-validation partitions data into $k$ equal-sized folds, trains on $k-1$ folds and validates on the remaining fold, repeating $k$ times. The average performance $\text{CV}_{(k)} = \frac{1}{k}\sum_{i=1}^{k} \text{MSE}_i$ provides a more robust estimate of generalization error. Stratified k-fold preserves class proportions in each fold, essential for imbalanced classification. Leave-one-out cross-validation (LOOCV) is the extreme case where $k=n$, providing near-unbiased but high-variance estimates.

K-Nearest Neighbors (KNN)

The K-Nearest Neighbors (KNN) algorithm is a non-parametric, instance-based learner that requires no explicit training phase. For a new query point $\mathbf{x}$, it identifies the $k$ training observations closest under a chosen distance metric, typically the Euclidean norm $d(\mathbf{x}, \mathbf{x}_i) = \|\mathbf{x} - \mathbf{x}_i\|_2$.

In classification, the predicted class $\hat{y}$ is the majority vote among the $k$ neighbors. In regression, the prediction is the (possibly inverse-distance-weighted) average of the neighbors' target values:

$$\hat{y}(\mathbf{x}) = \frac{\sum_{i \in \mathcal{N}_k(\mathbf{x})} \frac{y_i}{d(\mathbf{x},\mathbf{x}_i)}}{\sum_{i \in \mathcal{N}_k(\mathbf{x})} \frac{1}{d(\mathbf{x},\mathbf{x}_i)}}.$$

Common distance metrics include the Minkowski distance $d_p(x, x') = \left(\sum_{j=1}^m |x_j - x'_j|^p\right)^{1/p}$, which specializes to Manhattan ($p=1$) and Euclidean ($p=2$) distances. The hyper-parameter $k$ governs the bias–variance trade-off: small $k$ yields highly flexible, low-bias but high-variance boundaries; large $k$ smoothes the decision surface.

Although KNN is conceptually appealing, its performance degrades in high-dimensional spaces (the curse of dimensionality) because all points become roughly equidistant. Formally, as $m \to \infty$, $\frac{\max_i d(x_q, x_i) - \min_i d(x_q, x_i)}{\min_i d(x_q, x_i)} \xrightarrow{P} 0$, rendering the notion of "nearest" statistically meaningless. The algorithm's economic intuition is that of a pure market-follower: it predicts the outcome of a new observation by looking at the most similar historical cases, much like a real-estate appraisal that averages recent sale prices of comparable properties.

Decision Trees

A Decision Tree recursively partitions the input space $\mathcal{X}$ into axis-aligned hyper-rectangles, assigning a constant prediction (the majority class or the mean response) to each terminal leaf. The tree is grown by sequentially choosing the attribute and splitting threshold that maximize a measure of node purity.

Attribute Selection Criteria

For classification, given a set $S$ of samples and classes $c = 1,\dots,C$ with proportions $p_c$, the entropy captures the uncertainty:

$$H(S) = -\sum_{c=1}^{C} p_c \log_2 p_c.$$

The information gain (ID3 algorithm) of splitting on attribute $A$ is the reduction in entropy:

$$IG(S, A) = H(S) - \sum_{v \in \text{Values}(A)} \frac{|S_v|}{|S|} H(S_v).$$

Because information gain favours attributes with many distinct values, C4.5 introduces the gain ratio:

$$\text{GainRatio}(S,A) = \frac{IG(S,A)}{-\sum_{v} \frac{|S_v|}{|S|} \log_2 \frac{|S_v|}{|S|}}.$$

The CART algorithm, in contrast, uses the Gini impurity:

$$\text{Gini}(S) = 1 - \sum_{c=1}^{C} p_c^2,$$

which measures the probability of misclassification if a sample were labelled randomly according to the class distribution. For regression trees, splitting decisions are based on the reduction in mean squared error:

$$\text{SSE} = \sum_{i \in R_L} (y_i - \bar{y}_L)^2 + \sum_{i \in R_R} (y_i - \bar{y}_R)^2.$$

Tree-Building Algorithms and Pruning

• ID3 — restricted to categorical features, no missing-value handling, no pruning; it greedily maximizes information gain.
• C4.5 — extends ID3 to numeric attributes (candidate splits are midpoints between sorted values), handles missing values via fractional counts, and uses gain ratio. Post-pruning via error-based pruning simplifies the tree.
• CART — builds binary trees using Gini impurity (classification) or least-squares error (regression). It applies cost-complexity pruning: a sequence of subtrees $T_\alpha$ is generated by minimizing $R(T) + \alpha |T|$, where $R(T)$ is the misclassification/error rate and $|T|$ the number of leaves; $\alpha$ is selected by cross-validation.

From an economic perspective, decision trees mirror sequential investment decisions under uncertainty: each internal node represents a "test" that partitions the state space, and the leaf provides the expected payoff of the path. Their interpretability makes them attractive in regulated industries where model transparency is mandatory.

Random Forest and Ensemble Methods

Random Forest is an ensemble of decision trees trained via bagging (Bootstrap Aggregating). Each tree is trained on a bootstrap sample of the data, and at each split, only a random subset of $m$ features (typically $m \approx \sqrt{p}$ for classification, $p/3$ for regression) is considered. The final prediction is the majority vote (classification) or average (regression) across all trees. The random feature selection decorrelates the trees, reducing variance:

$$\text{Var}(\bar{f}) = \rho \sigma^2 + \frac{1-\rho}{B}\sigma^2,$$

where $\rho$ is the average pairwise correlation between trees, $\sigma^2$ is the variance of a single tree, and $B$ is the number of trees. Reducing $\rho$ through random feature subsets directly reduces ensemble variance.

Boosting trains models sequentially, with each subsequent model focusing on the errors of previous ones. The ensemble prediction is a weighted combination:

$$F(x) = \sum_{m=1}^{M} \alpha_m f_m(x).$$

AdaBoost adjusts sample weights: $w_i^{(m+1)} = w_i^{(m)} \exp(-\alpha_m y_i f_m(x_i))$, increasing the weight of misclassified points. Gradient Boosting frames boosting as functional gradient descent in function space. At each iteration $m$, a new tree $h_m(x)$ is fit to the negative gradient of the loss function:

$$r_{im} = -\left[\frac{\partial L(y_i, F(x_i))}{\partial F(x_i)}\right]_{F=F_{m-1}}.$$

XGBoost (Extreme Gradient Boosting) uses a second-order Taylor approximation of the loss:

$$\mathcal{L}^{(m)} \approx \sum_{i=1}^{n} \left[g_i f_m(x_i) + \frac{1}{2}h_i f_m(x_i)^2\right] + \Omega(f_m),$$

where $g_i = \partial_{\hat{y}^{(m-1)}} l(y_i, \hat{y}^{(m-1)})$ and $h_i = \partial^2_{\hat{y}^{(m-1)}} l(y_i, \hat{y}^{(m-1)})$ are the first and second order gradients. $\Omega(f) = \gamma T + \frac{1}{2}\lambda \sum_{j=1}^{T} w_j^2$ is a regularization term penalizing tree complexity. XGBoost includes built-in handling of missing values, column subsampling, and shrinkage (learning rate $\eta$).

LightGBM uses histogram-based algorithms and leaf-wise tree growth (expanding the leaf with the highest loss reduction) rather than level-wise growth. It introduces Gradient-based One-Side Sampling (GOSS) to retain instances with large gradients while randomly sampling those with small gradients, and Exclusive Feature Bundling (EFB) to bundle mutually exclusive features, dramatically reducing feature dimensionality.

Support Vector Machines (SVM)

Support Vector Machines (SVM) construct a maximum-margin separating hyperplane for binary classification. Given linearly separable data $\{(\mathbf{x}_i, y_i)\}_{i=1}^n$ with $y_i \in \{-1, +1\}$, the hard-margin SVM solves the convex quadratic program:

$$\min_{\mathbf{w}, b} \frac{1}{2}\|\mathbf{w}\|^2 \quad \text{s.t.} \quad y_i(\mathbf{w}^\top \mathbf{x}_i + b) \ge 1,\; i=1,\dots,n.$$

The Lagrangian dual reveals that the solution depends only on the subset of training points (the support vectors) with nonzero Lagrange multipliers $\alpha_i$, and the decision function becomes

$$f(\mathbf{x}) = \operatorname{sgn}\!\left(\sum_{i\in SV} \alpha_i y_i \mathbf{x}_i^\top \mathbf{x} + b\right).$$

When data are not perfectly separable, the soft-margin formulation introduces slack variables $\xi_i \ge 0$:

$$\min_{\mathbf{w}, b, \boldsymbol{\xi}} \frac{1}{2}\|\mathbf{w}\|^2 + C \sum_{i=1}^n \xi_i \quad \text{s.t.} \quad y_i(\mathbf{w}^\top \mathbf{x}_i + b) \ge 1 - \xi_i,\; \xi_i \ge 0.$$

The parameter $C$ controls the trade-off between margin width and classification error. The dual formulation is:

$$\max_{\alpha} \sum_{i=1}^n \alpha_i - \frac{1}{2} \sum_{i,j} \alpha_i \alpha_j y_i y_j \mathbf{x}_i^\top \mathbf{x}_j \quad \text{s.t.} \quad 0 \leq \alpha_i \leq C, \quad \sum_i \alpha_i y_i = 0.$$

Non-linear SVM and the Kernel Trick

Real-world data are rarely linearly separable. The kernel trick circumvents explicit feature mapping by replacing inner products with a symmetric, positive semi-definite kernel function $K(\mathbf{x}_i, \mathbf{x}_j) = \langle \phi(\mathbf{x}_i), \phi(\mathbf{x}_j) \rangle$. By Mercer's theorem, any continuous, symmetric, positive semi-definite kernel corresponds to an inner product in some reproducing kernel Hilbert space (RKHS). Common choices:

• Polynomial kernel: $K(\mathbf{x},\mathbf{z}) = (\gamma\,\mathbf{x}^\top \mathbf{z} + r)^d$
• Radial Basis Function (RBF) kernel: $K(\mathbf{x},\mathbf{z}) = \exp\!\left(-\gamma \|\mathbf{x} - \mathbf{z}\|^2\right)$
• Sigmoid kernel: $K(\mathbf{x},\mathbf{z}) = \tanh(\kappa \mathbf{x}^\top \mathbf{z} + c)$

The RBF kernel effectively lifts the data into an infinite-dimensional space, enabling highly flexible decision boundaries. The hyper-parameter $\gamma$ controls the influence of a single training example: large $\gamma$ leads to peaked Gaussians and complex, locally adapted boundaries prone to over-fitting.

Support Vector Regression (SVR)

Support Vector Regression (SVR) adapts the margin principle to regression. SVR fits a tube of radius $\epsilon$ around the predicted function $\hat{y} = \mathbf{w}^\top \mathbf{x} + b$ such that errors smaller than $\epsilon$ are ignored, while larger deviations incur a linear penalty:

$$\min_{\mathbf{w}, b, \boldsymbol{\xi}, \boldsymbol{\xi}^*} \frac{1}{2}\|\mathbf{w}\|^2 + C \sum_{i=1}^n (\xi_i + \xi_i^*)$$

subject to

$ \begin{aligned} y_i - \mathbf{w}^\top \mathbf{x}_i - b &\le \epsilon + \xi_i, \\ \mathbf{w}^\top \mathbf{x}_i + b - y_i &\le \epsilon + \xi_i^*, \\ \xi_i, \xi_i^* &\ge 0. \end{aligned} $

The $\epsilon$-insensitive loss function $L_\epsilon(y, \hat{y}) = \max(0, |y - \hat{y}| - \epsilon)$ creates a sparse solution where only points outside the $\epsilon$-tube contribute to the model.

Principal Component Analysis: Variance-Maximizing Projection

Principal Component Analysis (PCA) seeks an orthogonal basis along which the projected data exhibits maximal variance. Let $\mathbf{X} \in \mathbb{R}^{n \times d}$ be a mean-centered data matrix. The sample covariance matrix is

$$\mathbf{S} = \frac{1}{n-1}\mathbf{X}^\top \mathbf{X} \in \mathbb{R}^{d \times d}.$$

We wish to find a unit vector $\mathbf{w}_1$ that maximizes the variance of the projection $\mathbf{X}\mathbf{w}_1$:

$$\max_{\mathbf{w}} \; \mathbf{w}^\top \mathbf{S} \mathbf{w} \quad \text{subject to} \quad \mathbf{w}^\top \mathbf{w} = 1.$$

Forming the Lagrangian $\mathcal{L}(\mathbf{w}, \lambda) = \mathbf{w}^\top \mathbf{S}\mathbf{w} - \lambda(\mathbf{w}^\top\mathbf{w} - 1)$ and setting $\nabla_{\mathbf{w}} \mathcal{L} = \mathbf{0}$ yields the eigenvalue equation

$$\mathbf{S}\mathbf{w} = \lambda \mathbf{w}.$$

The variance-maximizing direction is the eigenvector associated with the largest eigenvalue $\lambda_1$. Subsequent components $\mathbf{w}_2, \mathbf{w}_3, \dots$ are obtained by repeating the optimization under the additional orthogonality constraint $\mathbf{w}_j^\top \mathbf{w}_i = 0$ for $j > i$, producing the ordered spectrum $\lambda_1 \geq \lambda_2 \geq \dots \geq \lambda_d \geq 0$.

Equivalently, PCA can be computed via the singular value decomposition (SVD) $\mathbf{X} = \mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^\top$, where the right singular vectors $\mathbf{V}$ are the principal directions and the singular values satisfy $\sigma_i^2 = (n-1)\lambda_i$. The $k$-dimensional representation is $\mathbf{Z} = \mathbf{X}\mathbf{V}_k = \mathbf{U}_k \boldsymbol{\Sigma}_k$, and the fraction of total variance retained is $\sum_{i=1}^k \lambda_i / \sum_{i=1}^d \lambda_i$ — a quantity routinely inspected via a scree plot to select $k$.

Practical implementation steps:

1. Center the data: subtract the mean of each feature.
2. (Optional) Scale to unit variance if features have different units.
3. Compute SVD of the centered (and scaled) matrix.
4. Project data onto the first $k$ right singular vectors.
5. Transform new data using the same centering/scaling and projection matrix $\mathbf{W}_k$.

from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

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

pca = PCA(n_components=0.95)  # retain 95% variance
X_pca = pca.fit_transform(X_scaled)

print(f"Explained variance ratio: {pca.explained_variance_ratio_}")
print(f"Number of components: {pca.n_components_}")

Intuitively, PCA identifies the axes of greatest spread in the data ellipsoid. In finance, the first principal component of a cross-section of asset returns often proxies the market factor; in macroeconomics, a handful of components summarize hundreds of indicators into "coincident" and "leading" indices. However, PCA has two structural limitations: it is sensitive to feature scaling (necessitating standardization when units differ), and it captures only second-order (linear correlation) structure. Two variables may be uncorrelated yet strongly dependent, a gap that ICA is designed to close.

Independent Component Analysis: Recovering Latent Sources

Independent Component Analysis (ICA) addresses a fundamentally different problem. Suppose we observe $d$ signals $\mathbf{x} = (x_1, \dots, x_d)^\top$ that are unknown linear mixtures of $d$ latent source signals $\mathbf{s} = (s_1, \dots, s_d)^\top$:

$$\mathbf{x} = \mathbf{A}\mathbf{s},$$

where $\mathbf{A} \in \mathbb{R}^{d \times d}$ is an unknown mixing matrix. The goal of ICA is to estimate an unmixing matrix $\mathbf{W} \approx \mathbf{A}^{-1}$ such that the recovered components $\mathbf{y} = \mathbf{W}\mathbf{x}$ are as statistically independent as possible. Two random variables $X$ and $Y$ are statistically independent if and only if $p(x, y) = p(x)\,p(y) \; \forall\, x, y$. This condition is strictly stronger than uncorrelatedness ($\text{Cov}(X,Y)=0$).

Why Non-Gaussianity Is the Key

A central result is that independence can only be identified when the sources are non-Gaussian. The argument rests on two pillars:

1. Rotational invariance of the Gaussian. If $\mathbf{s}$ is a vector of i.i.d. standard Gaussians, then any orthogonal rotation $\mathbf{Q}\mathbf{s}$ has the same distribution. The mixing matrix $\mathbf{A}$ is therefore unidentifiable from second-order statistics alone.
2. The Central Limit Theorem. A sum (or linear mixture) of independent non-Gaussian variables is more Gaussian than its constituents. Hence, searching for projections of $\mathbf{x}$ that are least Gaussian recovers the original sources.

This insight converts ICA into an optimization: find $\mathbf{w}$ such that $\mathbf{w}^\top \mathbf{x}$ is maximally non-Gaussian. Two standard measures are used:

• Kurtosis (the fourth standardized cumulant): $\kappa_4(y) = \mathbb{E}[y^4] - 3(\mathbb{E}[y^2])^2$. Gaussian variables have $\kappa_4 = 0$; super-Gaussian (sparse, heavy-tailed) sources have $\kappa_4 > 0$; sub-Gaussian (uniform-like) sources have $\kappa_4 < 0$.
• Negentropy, defined via the differential entropy $H(y) = -\int p(y)\log p(y)\,dy$: $J(y) = H(y_{\text{gauss}}) - H(y),$ where $y_{\text{gauss}}$ is a Gaussian with the same covariance as $y$. Since the Gaussian maximizes entropy for a given variance, $J(y) \geq 0$ with equality iff $y$ is Gaussian.

The FastICA Algorithm

The most widely used ICA solver, FastICA, maximizes an approximation to negentropy through a fixed-point iteration. After whitening $\mathbf{x}$ to $\mathbf{z}$ (so that $\mathbb{E}[\mathbf{z}\mathbf{z}^\top] = \mathbf{I}$, which PCA provides for free), one iterates for each component:

$$\mathbf{w}^+ = \mathbb{E}\!\left[\mathbf{z}\, g(\mathbf{w}^\top \mathbf{z})\right] - \mathbb{E}\!\left[g'(\mathbf{w}^\top \mathbf{z})\right]\mathbf{w}, \qquad \mathbf{w} \leftarrow \frac{\mathbf{w}^+}{\|\mathbf{w}^+\|},$$

where $g(u) = \tanh(u)$ or $g(u) = u\exp(-u^2/2)$ serves as a contrast function approximating negentropy. Convergence is cubic under mild conditions. Two intrinsic ambiguities remain: the order of the components and their sign are indeterminate; ICA recovers sources up to scaling and permutation.

from sklearn.decomposition import FastICA

ica = FastICA(n_components=3, random_state=42)
S_ica = ica.fit_transform(X_scaled)
A_ica = ica.mixing_  # estimated mixing matrix

Contrast with PCA: PCA yields uncorrelated components that maximize variance; ICA yields independent components that capture the data's non-Gaussian structure. In many signals, the interesting structure lives in the tails, not the bulk variance. PCA tends to capture high-variance Gaussian-like noise, while ICA can extract hidden factors, e.g., separating EEG artifacts from brain activity, or identifying latent factors in financial time series. For high-dimensional data, apply PCA to reduce to a manageable number of dimensions (e.g., retaining 99% variance) before ICA to improve stability and computational efficiency.