Statistics and Machine Learning Toolbox: Methods guide

Mathematical formulations, loss functions, optimization algorithms, cross-validation theory, diagnostic equations, and MATLAB capability comparison.

← Statistics and Machine Learning Toolbox overview · Getting started · Modeling guide · Examples · Reference

Mathematical formulations

The Statistics and Machine Learning Toolbox integrates classical statistical learning, regularized linear models, tree-based partitions, and unsupervised spectral decompositions. This guide provides exact mathematical specifications of all estimators, objectives, cross-validation algorithms, and diagnostic tests executed by the runtime.

---

1. L2-Regularized Logistic Regression

Logistic regression models the posterior probability of a binary response \(y \in \{0, 1\}\) conditioned on a feature vector \(\mathbf{x} \in \mathbb{R}^p\).

Logit link and probability formulation

The probability of the positive class (\(y = 1\)) is defined via the standard logistic sigmoid function \(\sigma(z)\):

\[P(y = 1 \mid \mathbf{x}; \mathbf{w}, b) = \sigma(\mathbf{w}^T \mathbf{x} + b) = \frac{1}{1 + e^{-(\mathbf{w}^T \mathbf{x} + b)}}\]

where \(\mathbf{w} \in \mathbb{R}^p\) is the weight vector and \(b \in \mathbb{R}\) is the intercept. The odds of the event occurring are:

\[\text{Odds}(\mathbf{x}) = \frac{P(y = 1 \mid \mathbf{x})}{1 - P(y = 1 \mid \mathbf{x})} = e^{\mathbf{w}^T \mathbf{x} + b}\]

Taking the natural logarithm yields the linear log-odds (logit):

\[\operatorname{logit}(P(y = 1 \mid \mathbf{x})) = \ln\left(\frac{P(y = 1 \mid \mathbf{x})}{1 - P(y = 1 \mid \mathbf{x})}\right) = \mathbf{w}^T \mathbf{x} + b\]

Penalized objective function

To prevent overfitting on collinear features, scikit-learn’s LogisticRegression minimizes the negative log-likelihood (binary cross-entropy) regularized by an L2 weight penalty:

\[\min_{\mathbf{w}, b} \; \frac{1}{2} \|\mathbf{w}\|_2^2 + C \sum_{i=1}^n \left[ - y_i \ln \sigma(\mathbf{w}^T \mathbf{x}_i + b) - (1 - y_i) \ln (1 - \sigma(\mathbf{w}^T \mathbf{x}_i + b)) \right]\]

where: - \(C > 0\) is the inverse regularization parameter specified on Model Config (Default = 1.0). Smaller values of \(C\) impose stronger regularization penalties. - \(\|\mathbf{w}\|_2^2 = \sum_{j=1}^p w_j^2\) is the squared Euclidean norm of the feature coefficients (the intercept \(b\) is unpenalized).

Optimization solver

The objective is strictly convex and twice-differentiable. Optimization is performed using the Limited-memory BFGS (L-BFGS) quasi-Newton algorithm. L-BFGS approximates the inverse Hessian matrix \(\mathbf{H}^{-1}\) using recent gradient updates without storing the full \(p \times p\) matrix, guaranteeing fast superlinear convergence in Pyodide.

Odds ratio interpretation

For standardized features (\(z_{i, j} = (x_{i, j} - \mu_j)/\sigma_j\)), the Odds Ratio (OR) associated with feature \(j\) represents the multiplicative change in the odds of churn per one standard deviation increase in that feature:

\[\text{OR}_j = \exp(w_j)\]

• If \(w_j > 0 \implies \text{OR}_j > 1\): Increasing feature \(j\) increases the likelihood of the event.
• If \(w_j < 0 \implies \text{OR}_j < 1\): Increasing feature \(j\) decreases the likelihood of the event.
• If \(w_j = 0 \implies \text{OR}_j = 1\): The feature has zero marginal predictive impact.

---

2. Ridge Regression (Tikhonov Regularization)

Ridge regression addresses ordinary least squares instability caused by multicollinear predictors by adding an L2 penalty to the residual sum of squares.

Objective function

Given target vector \(\mathbf{y} \in \mathbb{R}^n\) and standardized design matrix \(\mathbf{X} \in \mathbb{R}^{n \times p}\), Ridge regression solves:

\[\min_{\mathbf{w}} \; \frac{1}{2n} \|\mathbf{y} - \mathbf{X}\mathbf{w}\|_2^2 + \frac{\alpha}{2} \|\mathbf{w}\|_2^2\]

where \(\alpha \ge 0\) is the regularization weight configured on Model Config (Default = 1.0).

Analytical closed-form solution

Setting the gradient with respect to \(\mathbf{w}\) to zero:

\[\nabla_{\mathbf{w}} \mathcal{L} = -\frac{1}{n} \mathbf{X}^T (\mathbf{y} - \mathbf{X}\mathbf{w}) + \alpha \mathbf{w} = \mathbf{0}\]

yields the regularized normal equations:

\[\hat{\mathbf{w}} = (\mathbf{X}^T \mathbf{X} + n\alpha \mathbf{I}_p)^{-1} \mathbf{X}^T \mathbf{y}\]

SVD formulation and numerical stability

Inverting \(\mathbf{X}^T \mathbf{X}\) directly can be numerically ill-conditioned. Scikit-learn utilizes the Singular Value Decomposition (SVD) of \(\mathbf{X}\):

\[\mathbf{X} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^T\]

where \(\mathbf{\Sigma} = \operatorname{diag}(\sigma_1, \dots, \sigma_p)\) are the singular values. The Ridge estimator is computed stably as:

\[\hat{\mathbf{w}} = \mathbf{V} \operatorname{diag}\left( \frac{\sigma_j}{\sigma_j^2 + n\alpha} \right) \mathbf{U}^T \mathbf{y}\]

As \(\alpha \to 0\), the solution approaches the Moore-Penrose pseudo-inverse (OLS). As \(\alpha \to \infty\), weights shrink smoothly toward zero, eliminating coefficient explosion without forcing sparsity.

---

3. K-Means Clustering

K-Means is an unsupervised partition algorithm that segments \(n\) observations into \(K\) disjoint clusters \(\mathcal{S} = \{S_1, S_2, \dots, S_K\}\), where each cluster is parameterized by its mean centroid \(\boldsymbol{\mu}_k \in \mathbb{R}^p\).

Objective (Inertia)

K-Means minimizes the within-cluster sum of squared Euclidean distances (Inertia \(J\)):

\[\min_{\mathcal{S}, \{\boldsymbol{\mu}_k\}} \; J = \sum_{k=1}^K \sum_{\mathbf{x}_i \in S_k} \|\mathbf{x}_i - \boldsymbol{\mu}_k\|_2^2\]

Centroid seeding (k-means++)

To avoid suboptimal local minima, initial cluster centroids are chosen via the probabilistic k-means++ initialization algorithm: 1. Select the first centroid \(\boldsymbol{\mu}_1\) uniformly at random from the dataset. 2. For each sample \(\mathbf{x}_i\), compute the squared distance \(D(\mathbf{x}_i)^2\) to the nearest chosen centroid. 3. Choose the next centroid \(\boldsymbol{\mu}_k\) from the data points with probability proportional to \(D(\mathbf{x}_i)^2\): \[P(\mathbf{x}_i) = \frac{D(\mathbf{x}_i)^2}{\sum_{m=1}^n D(\mathbf{x}_m)^2}\] 4. Repeat steps 2 and 3 until \(K\) centroids have been initialized.

Lloyd’s algorithm iteration

The algorithm alternates between two steps until convergence: 1. Assignment step: Assign each sample \(\mathbf{x}_i\) to its nearest centroid: \[S_k^{(t)} = \left\{ \mathbf{x}_i : \|\mathbf{x}_i - \boldsymbol{\mu}_k^{(t)}\|_2^2 \le \|\mathbf{x}_i - \boldsymbol{\mu}_j^{(t)}\|_2^2 \; \forall j, 1 \le j \le K \right\}\] 2. Update step: Recalculate cluster centroids as the geometric mean of assigned points: \[\boldsymbol{\mu}_k^{(t+1)} = \frac{1}{|S_k^{(t)}|} \sum_{\mathbf{x}_i \in S_k^{(t)}} \mathbf{x}_i\]

Cluster validation metrics

• Silhouette Score (\(s\)): For sample \(i\), let \(a(i)\) be the mean intra-cluster distance to all other points in the same cluster, and let \(b(i)\) be the smallest mean distance to points in any other cluster. The sample silhouette is: \[s(i) = \frac{b(i) - a(i)}{\max(a(i), b(i))}, \quad s(i) \in [-1, 1]\] The global score is the average across all samples: \(\bar{s} = \frac{1}{n} \sum_{i=1}^n s(i)\). Values \(> 0.50\) indicate well-separated cluster geometry.
• Calinski-Harabasz Index (\(CH\)): Ratio of between-cluster dispersion to within-cluster dispersion: \[CH = \frac{\operatorname{Tr}(\mathbf{B}_K) / (K - 1)}{\operatorname{Tr}(\mathbf{W}_K) / (n - K)}\] where \(\mathbf{B}_K\) is the between-cluster scatter matrix and \(\mathbf{W}_K\) is the within-cluster scatter matrix. Higher values indicate tighter, more distinct clusters.

---

4. Decision Tree Classification (CART)

Decision trees construct a non-parametric classification model by recursively partitioning the feature space into hyper-rectangular regions.

Gini impurity

At node \(t\) containing \(N_t\) samples, the Gini impurity measures the expected error rate if labels were randomly drawn from the class distribution:

\[I_G(t) = 1 - \sum_{c=1}^C p(c \mid t)^2\]

where \(p(c \mid t)\) is the empirical proportion of class \(c\) in node \(t\), and \(C\) is the number of classes. For pure nodes, \(I_G(t) = 0\).

Splitting criterion and impurity reduction

For candidate feature \(j\) and split threshold \(\theta\), the data is split into left child \(t_L\) and right child \(t_R\). The split maximizes the impurity reduction (information gain):

\[\Delta I_G(t, j, \theta) = I_G(t) - \left( \frac{N_{t_L}}{N_t} I_G(t_L) + \frac{N_{t_R}}{N_t} I_G(t_R) \right)\]

Regularization and depth control

Unconstrained trees can expand until every leaf is pure (\(I_G = 0\)), memorizing training noise. The toolbox controls complexity through max_depth on Model Config (Default = 3), bounding maximum tree depth and ensuring robust out-of-fold generalization.

---

5. Principal Component Analysis (PCA)

PCA is an orthogonal linear transformation that maps \(p\)-dimensional data into a lower-dimensional coordinate system that maximizes feature variance.

Covariance matrix and spectral decomposition

For centered design matrix \(\mathbf{X}_c \in \mathbb{R}^{n \times p}\) (where \(\sum_i x_{i, j} = 0\)), the sample covariance matrix is:

\[\mathbf{\Sigma} = \frac{1}{n - 1} \mathbf{X}_c^T \mathbf{X}_c\]

Its spectral eigendecomposition is:

\[\mathbf{\Sigma} = \mathbf{V} \mathbf{\Lambda} \mathbf{V}^T\]

where \(\mathbf{V} = [\mathbf{v}_1, \dots, \mathbf{v}_p] \in \mathbb{R}^{p \times p}\) contains the orthogonal eigenvectors (principal component loadings), and \(\mathbf{\Lambda} = \operatorname{diag}(\lambda_1, \dots, \lambda_p)\) contains the non-negative eigenvalues arranged in descending order: \(\lambda_1 \ge \lambda_2 \ge \dots \ge \lambda_p \ge 0\).

Singular value decomposition (SVD)

To preserve numerical precision, scikit-learn computes PCA directly on \(\mathbf{X}_c\) using full SVD:

\[\mathbf{X}_c = \mathbf{U} \mathbf{S} \mathbf{V}^T\]

The singular values \(s_j\) relate to the covariance eigenvalues by:

\[\lambda_j = \frac{s_j^2}{n - 1}\]

Explained variance ratio

The proportion of total sample variance accounted for by the \(j\)-th principal component is:

\[\eta_j = \frac{\lambda_j}{\sum_{m=1}^p \lambda_m}\]

The cumulative explained variance for \(k\) components is \(\sum_{j=1}^k \eta_j\).

Coordinate projection and reconstruction

Observations are projected onto the top \(k\) principal component axes:

\[\mathbf{Z} = \mathbf{X}_c \mathbf{V}_k \in \mathbb{R}^{n \times k}\]

The reconstructed approximation in the original feature space is:

\[\hat{\mathbf{X}} = \mathbf{Z} \mathbf{V}_k^T = \mathbf{X}_c \mathbf{V}_k \mathbf{V}_k^T\]

The reconstruction error for sample \(i\) is \(\|\mathbf{x}_i - \hat{\mathbf{x}}_i\|_2\).

---

6. Stratified K-Fold Cross-Validation

Cross-validation provides an unbiased estimate of the generalization error of estimator \(f_\theta\).

flowchart TD
    D["Full Dataset D (n samples)"] --> Split["Stratified Partition into K Folds"]
    Split --> F1["Fold 1 (Val) | Folds 2..K (Train)"]
    Split --> F2["Fold 2 (Val) | Folds 1, 3..K (Train)"]
    Split --> FK["Fold K (Val) | Folds 1..K-1 (Train)"]
    F1 --> M1["Metric score M1"]
    F2 --> M2["Metric score M2"]
    FK --> MK["Metric score MK"]
    M1 --> Mean["CV Mean = (1/K) * sum(Mk)"]
    M2 --> Mean
    MK --> Mean

Mathematical formulation

The dataset \(D = \{(\mathbf{x}_i, y_i)\}_{i=1}^n\) is partitioned into \(K\) disjoint subsets \(D_1, D_2, \dots, D_K\) such that:

\[\bigcup_{k=1}^K D_k = D \quad \text{and} \quad D_j \cap D_k = \emptyset \; \forall j \ne k\]

In Stratified \(K\)-Fold, the partition preserves the class proportions in every fold:

\[\frac{|D_k \cap \{i : y_i = c\}|}{|D_k|} \approx \frac{|\{i : y_i = c\}|}{n} \quad \forall c \in \{0, 1\}, \forall k \in \{1, \dots, K\}\]

For fold \(k\), the model \(f_\theta^{(-k)}\) is fitted on the training subset \(D \setminus D_k\) and evaluated on the out-of-fold validation subset \(D_k\):

\[M_k = \text{Metric}(f_\theta^{(-k)}, D_k)\]

The cross-validation estimate is the arithmetic mean across all \(K\) folds:

\[\bar{M}_{\text{CV}} = \frac{1}{K} \sum_{k=1}^K M_k\]

---

7. The 7 Diagnostic Checks: Formal Definitions

The toolbox evaluates seven rigorous statistical health checks after every run:

Check 1: Data validation

Evaluates sample sufficiency and completeness: \[\text{Status} = \begin{cases} \text{PASS} & \text{if } n > p \text{ and } \sum_{i, j} \mathbf{1}_{\text{missing}}(x_{i, j}) = 0 \\ \text{WARN} & \text{if } n > p \text{ and } \sum_{i, j} \mathbf{1}_{\text{missing}}(x_{i, j}) > 0 \\ \text{FAIL} & \text{if } n \le p \end{cases}\]

Check 2: Convergence & configuration

Verifies solver termination status and hyperparameter boundaries: \[\text{Status} = \begin{cases} \text{PASS} & \text{if Status} = \text{"OPTIMAL"} \text{ and } \boldsymbol{\theta} \in \Omega_{\text{valid}} \\ \text{FAIL} & \text{otherwise} \end{cases}\]

Check 3: Model performance

Verifies that the primary metric satisfies domain viability thresholds: \[\text{Status}_{\text{classification}} = \begin{cases} \text{PASS} & \text{if Accuracy} \ge 0.80 \\ \text{WARN} & \text{if } 0.65 \le \text{Accuracy} < 0.80 \\ \text{FAIL} & \text{if Accuracy} < 0.65 \end{cases}\]

\[\text{Status}_{\text{regression}} = \begin{cases} \text{PASS} & \text{if } R^2 \ge 0.80 \\ \text{WARN} & \text{if } 0.60 \le R^2 < 0.80 \\ \text{FAIL} & \text{if } R^2 < 0.60 \end{cases}\]

\[\text{Status}_{\text{clustering}} = \begin{cases} \text{PASS} & \text{if Silhouette} \ge 0.40 \\ \text{WARN} & \text{if } 0.25 \le \text{Silhouette} < 0.40 \\ \text{FAIL} & \text{if Silhouette} < 0.25 \end{cases}\]

\[\text{Status}_{\text{PCA}} = \begin{cases} \text{PASS} & \text{if Cumulative ExpVar} \ge 0.80 \\ \text{WARN} & \text{if } 0.60 \le \text{Cumulative ExpVar} < 0.80 \\ \text{FAIL} & \text{if Cumulative ExpVar} < 0.60 \end{cases}\]

Check 4: Generalization & stability

Measures the absolute difference between training metric and cross-validation score: \[\Delta_{\text{gen}} = |M_{\text{train}} - \bar{M}_{\text{CV}}|\]

\[\text{Status} = \begin{cases} \text{PASS} & \text{if } \Delta_{\text{gen}} \le 0.15 \\ \text{WARN} & \text{if } 0.15 < \Delta_{\text{gen}} \le 0.25 \\ \text{FAIL} & \text{if } \Delta_{\text{gen}} > 0.25 \end{cases}\]

(For K-Means, checks that minimum cluster size \(\min |S_k| \ge 2\); for PCA, checks cumulative explained variance \(\ge 0.80\).)

Check 5: Multicollinearity & conditioning

Evaluates the 2-norm condition number of the design matrix \(\mathbf{X}\): \[\kappa(\mathbf{X}) = \frac{\sigma_{\max}(\mathbf{X})}{\sigma_{\min}(\mathbf{X})}\]

\[\text{Status} = \begin{cases} \text{PASS} & \text{if } \kappa(\mathbf{X}) < 30.0 \\ \text{WARN} & \text{if } 30.0 \le \kappa(\mathbf{X}) < 100.0 \\ \text{FAIL} & \text{if } \kappa(\mathbf{X}) \ge 100.0 \end{cases}\]

Check 6: Target balance & residual normality

• Classification: Evaluates minority class fraction \(\rho_{\min} = \frac{\min(n_0, n_1)}{n}\): \[\text{Status} = \begin{cases} \text{PASS} & \text{if } \rho_{\min} \ge 0.25 \\ \text{WARN} & \text{if } 0.10 \le \rho_{\min} < 0.25 \\ \text{FAIL} & \text{if } \rho_{\min} < 0.10 \end{cases}\]
• Regression: Evaluates residual normality via the Shapiro-Wilk test: \[W = \frac{\left( \sum_{i=1}^n a_i e_{(i)} \right)^2}{\sum_{i=1}^n (e_i - \bar{e})^2}\] \[\text{Status} = \begin{cases} \text{PASS} & \text{if } p\text{-value} \ge 0.05 \\ \text{WARN} & \text{if } p\text{-value} < 0.05 \end{cases}\]

Check 7: Outlier detection

Computes standardized residuals (or reconstruction errors for PCA / centroid distances for K-Means): \[z_i = \frac{r_i - \bar{r}}{\sigma_r}\]

Counts extreme observations where \(|z_i| > 2.5\): \[\text{Status} = \begin{cases} \text{PASS} & \text{if Count} = 0 \\ \text{WARN} & \text{if } 1 \le \text{Count} \le 2 \\ \text{FAIL} & \text{if Count} > 2 \end{cases}\]

---

MATLAB Statistics and Machine Learning Toolbox capability mapping

The table below provides a direct translation between MATLAB functions and the Boardflare toolbox implementation:

MATLAB Function / Workflow	MATLAB Capability	Boardflare Implementation	Notes & Equivalence
fitclinear(X, y, 'Learner', 'logistic')	Regularized binary classification	LogisticRegression(C=1.0, penalty='l2')	Scikit-learn L-BFGS solver with L2 weight shrinkage.
fitrlinear(X, y, 'Learner', 'leastsquares', 'Regularization', 'ridge')	Regularized linear regression	Ridge(alpha=1.0)	Exact closed-form Tikhonov regularization via SVD.
fitctree(X, y, 'MaxNumSplits', 2^depth)	CART classification tree	DecisionTreeClassifier(max_depth=3)	Gini impurity recursive partitioning.
kmeans(X, k, 'Distance', 'sqeuclidean')	K-Means clustering	KMeans(n_clusters=3, init='k-means++')	Identical k-means++ probabilistic initialization.
pca(X)	Principal component analysis	PCA(n_components=2)	Exact SVD of centered feature matrix \(\mathbf{X}_c\).
cvpartition(y, 'KFold', k)	Cross-validation fold generator	cross_val_score(cv=k)	Automated Stratified \(K\)-Fold and standard \(K\)-Fold.
predict(model, Xnew)	Model scoring	Predictions spill range (Predictions!A5:F55)	Real-time row-by-row prediction spill.
predict_proba / posterior scores	Posterior probability estimation	Predictions!D (Probability column)	Exact sigmoid class probabilities \(P(y=1)\).
silhouette(X, clust)	Cluster evaluation	Primary metric on Results + Check 3	Mean silhouette score across assigned clusters.
Classification Learner App	GUI model exploration	Start!B6 dropdown + Model Config	Interactive parameter tuning directly in the Excel grid.