Machine Learning with scikit-learn

Scikit-learn is the most widely used machine learning library in Python for classical (non-deep-learning) tasks. It provides a consistent, carefully designed API that covers the entire machine learning workflow: preprocessing raw data, training models, evaluating performance, tuning hyperparameters, and interpreting results. Engineers working with sensor data, quality control measurements, predictive maintenance logs, or simulation outputs will find scikit-learn's breadth of algorithms and utilities especially valuable. Understanding how the library is structured — and why it was designed that way — is the foundation for using it effectively.

The scikit-learn API Design

One of scikit-learn's greatest strengths is its deliberate uniformity. Regardless of whether you are training a simple linear model or a complex ensemble, the interface looks almost identical. This consistency is not accidental; it follows three core design principles: estimators expose a fit method, predictors expose a predict method, and transformers expose a transform method. Objects that both transform and predict (such as pipelines) combine these roles. This predictability dramatically reduces cognitive load and allows rapid experimentation.

The fit method is the training step. It accepts at minimum a two-dimensional feature array X (shape: n_samples × n_features) and, for supervised models, a one-dimensional target array y (shape: n_samples). Internally, fit adjusts the model's internal parameters — regression coefficients, tree split thresholds, support vectors, and so on — to best explain the relationship between X and y in the training data. After fit completes, the learned parameters are stored as attributes on the estimator object, typically indicated by a trailing underscore (for example, model.coef_ for a linear model). This trailing-underscore convention signals that an attribute was learned from data rather than set by the user.

from sklearn.linear_model import LinearRegression
import numpy as np

# Simulate temperature vs. pressure data
X_train = np.array([[100], [150], [200], [250], [300]])
y_train = np.array([1.2, 1.8, 2.5, 3.1, 3.9])

model = LinearRegression()
model.fit(X_train, y_train)

print(model.coef_)       # Learned slope, e.g. [0.0135]
print(model.intercept_)  # Learned intercept, e.g. -0.126

The predict method applies the fitted model to new input data to generate output predictions. It accepts the same kind of two-dimensional feature array as fit and returns a prediction array with one value per sample. Importantly, predict must only be called after fit; calling it before training raises a NotFittedError. For classifiers, predict returns the most likely class label. For regressors, it returns the estimated continuous value. Some classifiers also expose predict_proba, which returns the probability of each class rather than a hard label — useful when you need a calibrated confidence score alongside a decision.

# Predict pressure at new temperatures
X_new = np.array([[175], [225]])
predictions = model.predict(X_new)
print(predictions)  # e.g. [2.238, 2.913]

The score method provides a quick, default evaluation of the fitted model on a labeled dataset. For regressors, it returns the coefficient of determination R², where 1.0 is a perfect fit and values near or below 0 indicate the model performs no better than simply predicting the mean. For classifiers, it returns accuracy — the fraction of samples predicted correctly. The convenience of score makes it ideal for rapid sanity checks, but engineers should always consider whether the default metric is appropriate for their problem (for example, accuracy is misleading on imbalanced fault-detection datasets).

X_test = np.array([[120], [280]])
y_test = np.array([1.45, 3.65])
print(model.score(X_test, y_test))  # R² on held-out data

The profound benefit of this consistent interface is that swapping one algorithm for another requires almost no code changes — only the estimator class itself changes. This makes A/B comparisons between algorithms straightforward and reproducible:

from sklearn.ensemble import RandomForestRegressor
from sklearn.linear_model import Ridge

models = {
    "Linear":       LinearRegression(),
    "Ridge":        Ridge(alpha=1.0),
    "RandomForest": RandomForestRegressor(n_estimators=100, random_state=42),
}

for name, mdl in models.items():
    mdl.fit(X_train, y_train)
    print(f"{name}: R² = {mdl.score(X_test, y_test):.4f}")

Data Preprocessing and Pipelines

Raw engineering data almost never arrives in a form that machine learning models can consume directly. Sensor readings may span vastly different numerical ranges; log files may contain categorical identifiers; some measurements may be missing entirely. Preprocessing addresses these issues before training begins, and doing it correctly — especially in the context of model evaluation — is one of the most common sources of subtle bugs in machine learning projects.

StandardScaler transforms each feature so that it has zero mean and unit variance across the training set. This is achieved by subtracting the feature's mean and dividing by its standard deviation: z = (x − μ) / σ. Many algorithms are sensitive to the scale of input features. Gradient descent converges faster when features are on similar scales because the loss landscape becomes more spherical. Distance-based algorithms such as k-nearest neighbors and support vector machines produce results that depend directly on the Euclidean distances between samples — a feature measured in megapascals would otherwise dominate one measured in millimeters. Tree-based models are scale-invariant, but scaling rarely hurts them either.

from sklearn.preprocessing import StandardScaler
import numpy as np

# Two features: temperature (°C) and vibration amplitude (mm)
X = np.array([[80, 0.02], [100, 0.05], [120, 0.08], [140, 0.11]])

scaler = StandardScaler()
scaler.fit(X)           # Compute mean and std from training data
X_scaled = scaler.transform(X)

print(scaler.mean_)     # e.g. [110.0, 0.065]
print(scaler.scale_)    # e.g. [22.36, 0.033]
print(X_scaled)
# Each column now has mean ≈ 0 and std ≈ 1

A critical subtlety: the scaler must be fit only on training data and then applied (using transform, not fit_transform) to both training and test data. Fitting on the combined dataset causes data leakage — the model indirectly sees statistical properties of the test set during training, producing overly optimistic evaluation scores that do not generalize.

Categorical encoding is necessary because mathematical operations on raw category strings (such as "pump_A", "pump_B", "compressor") are undefined. The most common strategy is one-hot encoding, which creates a new binary column for each unique category. A feature with three categories becomes three binary columns, with exactly one column equal to 1 for each sample. Scikit-learn's OneHotEncoder handles this, and ColumnTransformer allows applying different transformers to different columns of a mixed dataset simultaneously.

from sklearn.preprocessing import OneHotEncoder
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import StandardScaler

# Mixed dataset: numerical column + categorical column
import pandas as pd
df = pd.DataFrame({
    "temperature": [80, 100, 120, 140],
    "component":   ["pump", "valve", "pump", "compressor"]
})

preprocessor = ColumnTransformer([
    ("num", StandardScaler(),    ["temperature"]),
    ("cat", OneHotEncoder(),     ["component"]),
])

X_processed = preprocessor.fit_transform(df)
print(X_processed.shape)  # (4, 4): 1 scaled + 3 one-hot columns

An important practical note: with k categories, one-hot encoding creates k columns, but some linear models require dropping one (the "dummy variable trap" or perfect multicollinearity). OneHotEncoder provides a drop='first' parameter for this purpose.

sklearn.pipeline.Pipeline solves a fundamental workflow problem: how to chain multiple preprocessing steps and the final model together so they behave as a single, coherent unit. A Pipeline accepts a list of (name, transformer) tuples, with the final element being the estimator. When you call pipeline.fit(X, y), it calls fit_transform on each transformer in sequence and then calls fit on the estimator using the transformed data. When you call pipeline.predict(X_new), it applies transform to each preprocessing step first and then passes the result to the estimator's predict method.

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import Ridge

pipe = Pipeline([
    ("scaler", StandardScaler()),
    ("model",  Ridge(alpha=1.0)),
])

pipe.fit(X_train, y_train)
predictions = pipe.predict(X_test)
score = pipe.score(X_test, y_test)

The most important benefit of using a pipeline during cross-validation is correctness. When scikit-learn's cross-validation tools (such as cross_val_score) split data into folds, a pipeline ensures that the scaler is refit on the training fold of each split, with the validation fold then transformed using that fold's statistics. Without a pipeline, if you scale the entire dataset before cross-validation, each fold's validation data has already influenced the scaler's parameters — a subtle but consequential form of leakage that inflates performance estimates.

from sklearn.model_selection import cross_val_score

scores = cross_val_score(pipe, X_all, y_all, cv=5, scoring="r2")
print(scores)        # R² for each of 5 folds
print(scores.mean()) # Mean cross-validated R²

Regression Models in scikit-learn

Regression predicts a continuous numeric output from input features. In engineering, typical regression targets include remaining useful life of a component, predicted yield of a process, energy consumption, or structural deflection under load. Scikit-learn provides a rich selection of regression algorithms that differ in assumptions, interpretability, computational cost, and predictive power.

LinearRegression fits a model of the form ŷ = β₀ + β₁x₁ + β₂x₂ + … + βₙxₙ by minimizing the sum of squared residuals (ordinary least squares, OLS). The solution is computed analytically (via matrix decomposition), making it exact and fast. Its primary virtue is interpretability: each coefficient βᵢ represents the expected change in the output per unit increase in feature xᵢ, holding all other features constant — a statement engineers can often explain to domain experts and present in technical reports. Its limitation is that it assumes a linear relationship between features and target, which may be a poor approximation for nonlinear physical systems.

from sklearn.linear_model import LinearRegression
model = LinearRegression()
model.fit(X_train, y_train)

# Inspect coefficients
for feature, coef in zip(feature_names, model.coef_):
    print(f"{feature}: {coef:.4f}")

Ridge and Lasso regression extend linear regression by adding a regularization penalty to the loss function, which discourages the model from fitting noise in the training data (overfitting). Ridge (L2 regularization) adds a penalty proportional to the sum of squared coefficients: loss = RSS + α Σ βᵢ². This shrinks all coefficients toward zero but rarely sets them exactly to zero, so all features remain in the model. Lasso (L1 regularization) uses a penalty proportional to the sum of absolute values of coefficients: loss = RSS + α Σ |βᵢ|. The L1 penalty has the geometric property of producing sparse solutions — many coefficients become exactly zero — effectively performing automatic feature selection. This is particularly useful in high-dimensional sensor datasets where only a subset of sensors truly predict the outcome.

from sklearn.linear_model import Ridge, Lasso

ridge = Ridge(alpha=10.0)
ridge.fit(X_train, y_train)
# All coefficients remain nonzero but shrunk

lasso = Lasso(alpha=0.1)
lasso.fit(X_train, y_train)
# Some coefficients exactly zero — those features eliminated
print((lasso.coef_ == 0).sum(), "features zeroed out")

The hyperparameter α controls the strength of regularization. Larger α increases the penalty and pushes coefficients further toward zero, increasing bias but reducing variance. The optimal α is typically found via cross-validated grid search (discussed later). RidgeCV and LassoCV are convenient variants that perform this search internally.

Ensemble regressors such as RandomForestRegressor aggregate predictions from many individual decision trees, each trained on a random bootstrap sample of the data and considering only a random subset of features at each split. This combination of bagging (bootstrap aggregating) and feature randomization reduces variance substantially relative to a single decision tree, while retaining the ability to model complex nonlinear relationships. Random forests require little preprocessing (no scaling needed), handle mixed feature types gracefully, and are robust to outliers. Their main downsides are computational cost (especially for large numbers of trees) and reduced interpretability compared to linear models.

from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor

rf = RandomForestRegressor(n_estimators=200, max_depth=10, random_state=42)
rf.fit(X_train, y_train)
print("RF R²:", rf.score(X_test, y_test))

gb = GradientBoostingRegressor(n_estimators=200, learning_rate=0.05, random_state=42)
gb.fit(X_train, y_train)
print("GB R²:", gb.score(X_test, y_test))

The principle of comparing multiple models on a held-out test set before making a final selection is fundamental to sound engineering practice. The held-out test set — data that was never seen during training or hyperparameter tuning — provides an unbiased estimate of how the model will perform on genuinely new data in production. Evaluating on training data alone will almost always produce optimistic and misleading results due to overfitting. A standard workflow partitions data into training (60–70%), validation or cross-validation (15–20%), and test (15–20%) portions, with the test set locked away until the final evaluation.

Classification Models in scikit-learn

Classification predicts a discrete class label. Engineering applications include fault detection (normal vs. faulty), quality classification (pass/fail/rework), anomaly detection in process control, and categorization of operational modes. Scikit-learn's classifiers share the same API as regressors, making it straightforward to benchmark several approaches.

LogisticRegression is, despite its name, a classification algorithm. It models the log-odds of class membership as a linear combination of features and passes the result through the logistic (sigmoid) function to produce a probability between 0 and 1. For binary classification: P(y=1|x) = 1 / (1 + e^(−(β₀ + β₁x₁ + … + βₙxₙ))). A threshold (default 0.5) converts this probability to a class label. Logistic regression is highly interpretable — its coefficients indicate how each feature shifts the log-odds of belonging to the positive class — and it is fast to train. It can also be extended to multi-class problems using one-vs-rest or softmax strategies. Like linear regression, it benefits from feature scaling and supports L1/L2 regularization via the penalty parameter.

from sklearn.linear_model import LogisticRegression
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split

X, y = make_classification(n_samples=500, n_features=10, random_state=42)
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.2, random_state=42)

clf = LogisticRegression(C=1.0, max_iter=1000)
clf.fit(X_tr, y_tr)
print("Accuracy:", clf.score(X_te, y_te))
print("Probabilities:", clf.predict_proba(X_te[:3]))

Decision tree classifiers partition the feature space into rectangular regions using a series of binary splits chosen to maximize a purity criterion such as Gini impurity or information gain. The resulting model can be visualized as a flowchart of if-then rules, which is an especially valuable property in engineering contexts: a maintenance engineer can follow the tree's logic manually, verify that it aligns with domain knowledge, and explain decisions to non-technical stakeholders. Scikit-learn's export_graphviz or plot_tree functions render the tree visually. The main risk is overfitting — an unconstrained tree will perfectly memorize training data. Controlling max_depth, min_samples_split, and min_samples_leaf prevents this.

from sklearn.tree import DecisionTreeClassifier, plot_tree
import matplotlib.pyplot as plt

dt = DecisionTreeClassifier(max_depth=4, random_state=42)
dt.fit(X_tr, y_tr)
print("DT Accuracy:", dt.score(X_te, y_te))

# Visualize (requires matplotlib)
fig, ax = plt.subplots(figsize=(20, 8))
plot_tree(dt, feature_names=[f"feat_{i}" for i in range(10)],
          class_names=["normal","fault"], filled=True, ax=ax)
plt.savefig("decision_tree.png", dpi=150, bbox_inches="tight")

Support Vector Machines (SVMs) find the hyperplane in feature space that maximizes the margin — the distance between the hyperplane and the nearest training samples from each class (the support vectors). This maximum-margin principle provides good generalization even in high-dimensional spaces, which is why SVMs were historically popular for text classification and genomics where features far outnumber samples. The kernel trick allows SVMs to implicitly operate in very high-dimensional or even infinite-dimensional feature spaces without explicitly computing coordinates there. Common kernels include the radial basis function (RBF), polynomial, and sigmoid kernels. SVMs are sensitive to feature scale, so StandardScaler is nearly always necessary. Their main disadvantages are computational cost on large datasets and difficulty interpreting the final model.

from sklearn.svm import SVC
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler

svm_pipe = Pipeline([
    ("scaler", StandardScaler()),
    ("svm",    SVC(kernel="rbf", C=10.0, gamma=0.1, probability=True)),
])
svm_pipe.fit(X_tr, y_tr)
print("SVM Accuracy:", svm_pipe.score(X_te, y_te))

Class imbalance is a pervasive problem in engineering fault-detection datasets. In real facilities, faults are rare events; a dataset might contain 950 normal samples and only 50 fault samples. A naive classifier that always predicts "normal" achieves 95% accuracy — a misleadingly high score that reflects no useful learning. Scikit-learn provides several remedies. The class_weight='balanced' parameter (available in LogisticRegression, SVC, DecisionTreeClassifier, and others) automatically adjusts the contribution of each sample to the loss function inversely proportional to class frequency, effectively penalizing misclassification of the minority class more heavily. The imbalanced-learn library (compatible with scikit-learn) provides oversampling algorithms such as SMOTE (Synthetic Minority Oversampling Technique), which generates synthetic minority samples by interpolating between existing ones. Undersampling removes majority class samples. In imbalanced settings, evaluation should use metrics such as F1-score, precision-recall AUC, or the Matthews Correlation Coefficient rather than accuracy.

from sklearn.linear_model import LogisticRegression
from sklearn.metrics import classification_report

# Use balanced class weights for imbalanced fault detection
clf_balanced = LogisticRegression(class_weight='balanced', max_iter=1000)
clf_balanced.fit(X_tr, y_tr)
y_pred = clf_balanced.predict(X_te)

# Detailed per-class metrics
print(classification_report(y_te, y_pred, target_names=["normal", "fault"]))

Hyperparameter Tuning

Every machine learning algorithm exposes hyperparameters — configuration choices that govern how the algorithm trains but are not learned from data themselves. Examples include the regularization strength α in Ridge regression, the number of trees in a random forest, the kernel type and margin parameter C in an SVM, and the learning rate in gradient boosting. Choosing good hyperparameters is crucial: poor choices can mean the difference between a model that generalizes well and one that dramatically overfits or underfits. Since the right hyperparameters depend on the specific dataset, they must be selected empirically.

Cross-validated grid search (GridSearchCV) is the systematic approach. You define a parameter grid — a dictionary mapping hyperparameter names to lists of candidate values — and GridSearchCV trains and evaluates every possible combination using k-fold cross-validation. The combination that produces the best mean cross-validated score is selected. After the search completes, the best model is retrained on the full training set. GridSearchCV wraps any estimator or pipeline and exposes the same fit/predict/score API, so it can itself be used as an estimator.

from sklearn.model_selection import GridSearchCV
from sklearn.ensemble import RandomForestClassifier

param_grid = {
    "n_estimators": [50, 100, 200],
    "max_depth":    [None, 5, 10, 20],
    "min_samples_split": [2, 5, 10],
}

rf = RandomForestClassifier(random_state=42)
grid_search = GridSearchCV(
    estimator=rf,
    param_grid=param_grid,
    cv=5,
    scoring="f1",
    n_jobs=-1,       # Use all CPU cores
    verbose=1,
)
grid_search.fit(X_tr, y_tr)

print("Best params:", grid_search.best_params_)
print("Best CV F1:", grid_search.best_score_)
print("Test F1:", grid_search.score(X_te, y_te))

The parameter grid above has 3 × 4 × 3 = 36 combinations. With 5-fold cross-validation, this means 180 model fits. For large datasets or models with many parameters, this becomes computationally prohibitive — which is where randomized search offers a practical advantage.

RandomizedSearchCV samples a specified number of hyperparameter combinations at random from the defined parameter distributions rather than exhaustively evaluating all combinations. A key insight from research (Bergstra & Bengio, 2012) is that random search is more efficient than grid search when only a few hyperparameters strongly affect performance — grid search wastes evaluations on redundant combinations along irrelevant dimensions, whereas random search spreads exploration more broadly. Distributions can be specified using scipy probability distributions (e.g., randint, uniform), allowing continuous sampling rather than discrete lists.

from sklearn.model_selection import RandomizedSearchCV
from scipy.stats import randint, uniform

param_dist = {
    "n_estimators":      randint(50, 500),
    "max_depth":         [None, 5, 10, 20, 30],
    "min_samples_split": randint(2, 20),
    "max_features":      uniform(0.3, 0.7),  # Fraction of features
}

rand_search = RandomizedSearchCV(
    estimator=RandomForestClassifier(random_state=42),
    param_distributions=param_dist,
    n_iter=50,         # Evaluate 50 random combinations
    cv=5,
    scoring="f1",
    n_jobs=-1,
    random_state=42,
)
rand_search.fit(X_tr, y_tr)
print("Best params:", rand_search.best_params_)

The principle that hyperparameter tuning must never use the test set cannot be overstated. If you evaluate many hyperparameter combinations against the test set and pick the best, you are effectively fitting the hyperparameters to the test set, causing optimistically biased performance estimates. The test set must remain completely unseen until the very final evaluation of the single selected model. Cross-validation on the training set provides an unbiased proxy for generalization performance during the tuning phase.

Good engineering practice also demands documenting the tuning process. This includes recording the hyperparameter search space explored, the scoring metric used, the cross-validation strategy, the winning combination, and the final test set performance. This documentation supports reproducibility (a colleague or auditor can re-run the search and obtain the same result given the same random_state seeds), traceability (demonstrating that a safety-critical model was selected through a rigorous, systematic process), and future refinement (knowing what was already tried prevents redundant work when revisiting the project).

Model Interpretability and Feature Importance

Predictive accuracy alone is often insufficient in engineering practice. A regulatory agency certifying a predictive maintenance system, a process engineer validating a quality model, or a safety manager approving an anomaly detection system all need to understand why a model makes its predictions — not just how accurate it is. Interpretability techniques bridge the gap between model outputs and domain understanding.

Feature importance from tree-based models is the most accessible interpretability tool. During the training of a decision tree or random forest, every split in every tree is associated with a reduction in impurity (Gini impurity for classifiers, mean squared error for regressors). Summing these impurity reductions across all splits that involve a given feature, weighted by the number of samples passing through each node and normalized across all features, gives that feature's importance score. In a random forest, these values are averaged across all trees. The result is a score between 0 and 1 for each feature, with higher values indicating greater contribution to the model's predictions.

import pandas as pd
import matplotlib.pyplot as plt

feature_names = [f"sensor_{i}" for i in range(10)]
importances = rf.feature_importances_

importance_df = pd.DataFrame({
    "feature":    feature_names,
    "importance": importances,
}).sort_values("importance", ascending=False)

print(importance_df)

importance_df.plot(kind="bar", x="feature", y="importance", legend=False)
plt.ylabel("Gini Importance")
plt.title("Feature Importances — Random Forest")
plt.tight_layout()
plt.savefig("feature_importances.png")

A known limitation of impurity-based importances is their bias toward high-cardinality features (features with many unique values), since such features offer more split opportunities and can appear important even if their relationship with the target is weak. This is where permutation importance provides a more reliable alternative.

Permutation importance is model-agnostic: it can be applied to any fitted model, not just tree-based ones. The procedure is: (1) evaluate the model's performance on a validation set; (2) randomly shuffle the values of one feature across all samples, breaking its relationship with the target while leaving all other features intact; (3) re-evaluate performance; (4) the drop in performance is the permutation importance of that feature. A large drop means the model relied heavily on that feature. A negligible drop means the feature contributed little (or that another correlated feature compensated). This is repeated for each feature, and optionally repeated multiple times with different shuffles to produce a confidence interval.

from sklearn.inspection import permutation_importance

result = permutation_importance(
    rf,
    X_te, y_te,
    n_repeats=30,
    random_state=42,
    scoring="f1",
)

perm_df = pd.DataFrame({
    "feature": feature_names,
    "importance_mean": result.importances_mean,
    "importance_std":  result.importances_std,
}).sort_values("importance_mean", ascending=False)

print(perm_df)

If impurity-based and permutation importance give substantially different rankings, it is often worth investigating whether correlated features are masking each other's contributions — a common situation when multiple sensors measure related physical quantities.

Partial dependence plots (PDPs) visualize the marginal effect of one or two features on the model's predicted outcome, averaging out the effects of all other features. For a chosen feature xⱼ, the PDP is computed by: (1) for each unique value in a grid spanning xⱼ's range, replace all samples' values of xⱼ with that grid value; (2) run the modified dataset through the model; (3) average the predictions. The resulting curve shows how the predicted output changes as xⱼ varies, marginalizing over the distribution of all other features. This can reveal whether a relationship is linear, monotonic, U-shaped, or has threshold effects — physical insights that complement the quantitative accuracy metrics.

from sklearn.inspection import PartialDependenceDisplay

fig, ax = plt.subplots(figsize=(12, 4))
PartialDependenceDisplay.from_estimator(
    rf,
    X_tr,
    features=[0, 2, (0, 2)],   # Feature 0, feature 2, and their interaction