Building models with scikit-learn

scikit-learn is where machine learning on a table of numbers happens in Python. Linear models, trees, forests, boosting, support-vector machines, clustering, dimensionality reduction, and around them the parts that matter more than the algorithm: preprocessing, cross-validation, hyperparameter search, metrics. It began as a Google Summer of Code project in 2007 and became what it is when a group at INRIA in France took it over in 2010 (about). The Python stack post placed it among its neighbours; this post is about using it well.

The example is the one the ML basics series used: five tickers, daily features computed from prices (recent returns, recent volatility, distance from the moving average), and the question will the next 20 trading days be a stress month, volatility above 30% a year. Train on 2019 to 2021, test on 2022 to 2023. Every number below comes from one script on that data.

The one idea: everything is an estimator

scikit-learn has one interface, and once you know it you know the library. Every object is an estimator: you construct it with its settings (hyperparameters), you call fit(X, y) to learn from data, and afterwards it has learned attributes whose names end in an underscore (coef_, feature_importances_, classes_). Beyond that, three kinds:

The design was a deliberate choice, written up by the core developers in 2013 (Buitinck et al.): a shared interface for every unit, composition over inheritance, and sensible defaults so that a model works before it is tuned. Géron’s Hands-On Machine Learning summarises the same principles as consistency, inspection, non-proliferation of classes, composition and sensible defaults. The practical consequence is that swapping a model is a one-line change, and the code around it (splitting, scaling, scoring, searching) does not care which model it holds.

from sklearn.compose import ColumnTransformer
from sklearn.dummy import DummyClassifier
from sklearn.ensemble import HistGradientBoostingClassifier, RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import OneHotEncoder, StandardScaler

FEATURES = ["ret_1d", "ret_5d", "ret_20d", "vol_5d", "vol_20d", "sma_ratio", "bollinger"]

def prep():                       # a fresh preprocessor per pipeline (see below)
    return ColumnTransformer([("num", StandardScaler(), FEATURES),
                              ("ticker", OneHotEncoder(handle_unknown="ignore"), ["ticker"])])

models = {
    "always 'calm'": DummyClassifier(strategy="prior"),
    "logistic regression": Pipeline([
        ("prep", prep()), ("model", LogisticRegression(max_iter=1000))]),
    "random forest": Pipeline([
        ("prep", prep()), ("model", RandomForestClassifier(
            n_estimators=300, max_features=3, min_samples_leaf=20, random_state=0))]),
    "gradient boosting": Pipeline([
        ("prep", prep()), ("model", HistGradientBoostingClassifier(
            max_depth=3, learning_rate=0.05, max_iter=200, random_state=0))]),
}
for name, m in models.items():
    m.fit(X_train, y_train)
    p = m.predict_proba(X_test)[:, 1]
    print(name, m.score(X_test, y_test), roc_auc_score(y_test, p))

Four models, one loop. Two things in that table before going further. First, the dummy classifier that always answers “calm” scores 77% accuracy, because 77% of test months are calm; accuracy is the wrong score for this problem, and the next section says what to use. Second, the plain logistic regression is the best of the three real models on this data, and by a margin that holds (the boosting post found the same on a sister problem). A linear model with scaled inputs is the baseline to beat. Do not skip it.

Preprocessing belongs inside the pipeline

The ColumnTransformer above scales the seven numeric columns and one-hot encodes the ticker, and the Pipeline glues that to the model. The glue is not there for tidiness. A transformer fitted on the whole dataset has seen the test set. StandardScaler learns a mean and a standard deviation; fit it on all the rows and the test rows’ statistics leak into training. Inside a pipeline, fit runs the scaler on the training rows only and predict reuses those numbers on the test rows, and the same holds inside every cross-validation fold. The same reasoning covers SimpleImputer (the fill value for missing data is learned) and OneHotEncoder (the set of categories is learned; handle_unknown="ignore" is what keeps a ticker that was not in training from crashing prediction).

Two things the library does not do for you:

• Estimators are stateful, and a pipeline fits its steps in place. The first version of the script built one ColumnTransformer and passed the same object to four pipelines. Each fit refitted it, which was harmless while every pipeline saw the same columns, and then a fifth pipeline fed it rescaled columns and silently broke the other four: their test AUC fell from $\key{0.83}$ to $\key{0.77}$ with no error. Build a fresh object per pipeline, or use sklearn.base.clone. Meta-estimators such as GridSearchCV and cross_val_score clone for you; Pipeline does not.
• Custom steps are a small class. Anything with fit and transform can go in a pipeline; inheriting from BaseEstimator and TransformerMixin gives it fit_transform, get_params and a sensible repr for free. Here is one that turns daily volatility into the yearly percentage the rest of the site quotes, so that the numbers mean something when the model is inspected:

from sklearn.base import BaseEstimator, TransformerMixin

class Annualise(BaseEstimator, TransformerMixin):
    def __init__(self, columns=("vol_5d", "vol_20d")):
        self.columns = columns        # stored as given, never changed

    def fit(self, X, y=None):
        return self                   # nothing to learn

    def transform(self, X):
        X = X.copy()
        X[list(self.columns)] = X[list(self.columns)] * (252 ** 0.5) * 100
        return X

pipe = make_pipeline(Annualise(), prep(), LogisticRegression(max_iter=1000))
pipe.fit(X_train, y_train)

The constructor stores its arguments untouched; that is the one rule the library’s developer guide insists on, because get_params and set_params (which grid search uses) read and write those attributes by name. The model with this step scores the same $\key{0.831}$ as before, since the scaler undoes a constant factor; the step is there so that vol_20d reads as “31% a year” rather than “0.0195” when you look inside.

Scores: accuracy is the wrong one here

The confusion matrix of the logistic regression on the test months, at the default cut-off of 0.5:

TP, stress months caughtFN, stress months missedFP, false alarms

$$\begin{aligned} \text{precision} &= \frac{\good{TP}}{\good{TP} + \risk{FP}} = \frac{132}{132 + 80} = 0.62 \\[10pt] \text{recall} &= \frac{\good{TP}}{\good{TP} + \bad{FN}} = \frac{132}{132 + 412} = 0.24 \end{aligned}$$

When the model says “stress” it is right $\key{62\%}$ of the time, and it catches only $\key{24\%}$ of the stress months. The 0.5 cut-off is a default, not a decision. A risk desk that wants to catch most stress months accepts more false alarms and moves the cut-off down. class_weight="balanced" does that at training time by weighting the rare class up: recall rises to $\key{70\%}$ (381 of 544 caught) and precision falls to $\key{45\%}$, with the AUC barely moved (0.831 to 0.829), because AUC measures the ranking of the probabilities, not the cut-off. That is why the tables here report AUC and average precision instead of accuracy: both score how well the model separates the classes, whatever cut-off is chosen later. classification_report prints precision, recall and F1 per class in one call; use it before quoting an accuracy.

Cross-validation: do not shuffle time

Cross-validation is how a score is estimated without touching the test set: split the training rows into folds, fit on all but one, score on the held-out one, repeat. The default splitter, KFold with shuffling, assigns rows to folds at random. On a time series that is a leak. Rows from the same week share almost everything (the same market, nearly the same 20-day volatility), so a model validated on a random Tuesday has trained on the Monday and Wednesday around it. The honest splitter is TimeSeriesSplit: every fold trains on the past and validates on the block that follows.

from sklearn.model_selection import KFold, TimeSeriesSplit, cross_val_score

gb = models["gradient boosting"]
shuffled = KFold(5, shuffle=True, random_state=0)
cross_val_score(gb, X_train, y_train, cv=shuffled, scoring="roc_auc")
cross_val_score(gb, X_train, y_train, cv=TimeSeriesSplit(5), scoring="roc_auc")

Panel 1 of the chart shows the two. The shuffled folds agree with each other and promise $\key{0.91}$; the model then delivers $\key{0.82}$ on real future data. The time-ordered folds are noisier and less flattering, and that is useful information: fold 2 trained on a calm 2019 and was asked about early 2020, and scored $\key{0.43}$, worse than a coin flip. A model that has never seen a crash does not know what one looks like. Shuffling hides exactly that.

Note. The rule is older than scikit-learn: validate on data that comes after the data you trained on, the way the real test will. The library will happily shuffle time series for you, because KFold has no idea the rows are dated. Choosing the splitter is your job, and it is the single setting that most changes what a cross-validated score means.

Grid search, on the right splitter

GridSearchCV takes an estimator, a dictionary of settings to try, a splitter and a scoring rule, fits every combination on every fold, and keeps the best. The model__learning_rate syntax (double underscore) reaches into a pipeline step by name, which is why pipelines and grid search fit together. The search below tries 12 combinations for the boosting model, and a second search tries 11 values of the regularisation strength C for the logistic regression.

from sklearn.model_selection import GridSearchCV

settings = {"model__learning_rate": [0.01, 0.03, 0.1, 0.3],
            "model__max_depth": [2, 3, 5]}
grid = GridSearchCV(gb, settings, cv=TimeSeriesSplit(5), scoring="roc_auc", n_jobs=-1)
grid.fit(X_train, y_train)
grid.best_params_, grid.best_score_       # chosen on validation folds only
p = grid.best_estimator_.predict_proba(X_test)[:, 1]
roc_auc_score(y_test, p)                  # looked at once, after the choice

Panel 2 shows the boosting grid: the shallowest, slowest-learning trees won, which is the usual answer on small noisy data, and the whole surface spans 0.73 to 0.79 with fold standard deviations around 0.1, so the differences between neighbours are inside the noise. The tuned boosting model reaches $\key{0.822}$ on the test set, up from 0.819 untuned; the tuned logistic regression reaches $\key{0.832}$. The test numbers were looked at once, after the choice was made. That is what makes them test numbers. cv_results_ holds every combination’s mean and spread, and reading it is more useful than reading best_params_.

Two habits from the same section of the library. RandomizedSearchCV samples settings instead of trying all of them, and finds as good an answer in a fraction of the fits once the grid has more than two or three dimensions. And n_jobs=-1 runs the folds on every core; scikit-learn has no GPU support, and parallel folds is where its speed comes from.

Inspection: what did it learn?

The learned attributes with the trailing underscore are the first place to look (coef_ for a linear model, feature_importances_ for trees). The better tool is permutation importance: shuffle one column of the test set, score again, and see how much was lost. For the tuned boosting model, shuffling the 20-day volatility costs $\risk{0.215}$ of AUC; no other column costs more than $\risk{0.004}$. The model is close to a function of one feature: this month’s volatility predicts next month’s. The random forests post got too. Panel 3 of the chart shows the three models’ ROC curves on the test months, close enough that the choice between them comes down to simplicity, and the simplest wins.

Saving the model

pickle and joblib write any fitted estimator to disk, and the library’s persistence guide gives the two warnings that matter. There is no supported way to load a model into a different version of scikit-learn, so record the version next to the file. And a pickle executes code when loaded, so never load one from a source you do not trust; skops is the format for that case, and ONNX for serving predictions without Python. Save the training script, the data snapshot and the cross-validation score with the model. Without them the file is a number with no history.

Where it stops

scikit-learn works on a table that fits in memory, on CPU, with classical methods. That covers more real problems than the attention paid to the alternatives suggests. The Python stack post says where to go when it is not enough: XGBoost and LightGBM for the last few points of accuracy on a table (both plug into pipelines and grid search unchanged), PyTorch when the input is images, text or audio, Dask or Spark when the table does not fit one machine. Two more limits are not in the manual: there is no time-series forecasting in the library (you build the lagged features yourself, as this example did, and then it is a table), and the defaults are tuned for independent rows, so every splitter, every shuffle and every importance measure needs a second look when the rows are dated.

Why banks use it

• The interface is auditable: a pipeline is a list of named steps with stored settings, and get_params() prints the whole thing. A model-risk reviewer can read it.
• Determinism: random_state everywhere, so a result reproduces.
• Linear models and trees, the methods a regulator will accept an explanation for, are first-class, with the inspection tools beside them.
• It runs on the machine the data is already on. No cluster, no GPU, no new infrastructure to approve.

Summary

Next

The ML basics series walks through each model this post used, from linear regression to boosting, on the same data. The backtesting post covers the same leakage and look-ahead traps in their trading form.

References

• Aurélien Géron, Hands-On Machine Learning with Scikit-Learn, Keras, and TensorFlow (2nd ed., O’Reilly, 2019), ch. 2 (the estimator design, transformation pipelines, custom transformers, grid search) and ch. 3 (confusion matrix, precision and recall, ROC curves). The author’s companion notebooks are at github.com/ageron/handson-mlp.
• James, Witten, Hastie, Tibshirani & Taylor, An Introduction to Statistical Learning, with Applications in Python (Springer, 2023), 5.1 on cross-validation.
• Lars Buitinck et al., “API design for machine learning software: experiences from the scikit-learn project” (2013), arXiv:1309.0238.
• scikit-learn user guide, the four pages this post leans on: pipelines and composite estimators, cross-validation, developing estimators and model persistence.