Every image and example was taken from the following sources:

• Lecture 7 — Seleção de Modelos Lineares e Regularização — Machine Learning 2026 Summer Course at IMPA
• Machine Learning 2026, Prof. Paulo Orenstein

Introduction

These notes cover Lecture 7 of the IMPA Summer Course — Machine Learning 2026 (Paulo Orenstein), one article per lecture. You can find the articles for the earlier lectures on my Medium profile.

Everything below is self-contained: you do not need to read the other posts to follow this one, as long as you are comfortable with least squares, the bias-variance trade-off, cross-validation and a little linear algebra.

The lecture covers Chapter 6 of the course textbook and asks what to do when a linear model has too many features. It starts with three open questions, then works through three families of answers: subset selection (best subset, forward and backward stepwise), shrinkage (ridge and the lasso, including how they relate to best subset and when each one wins) and dimension reduction (principal components regression).

Three open questions, one phenomenon

Linear models are powerful: interpretable, extensible and good at prediction. The idea is simple, but understanding how each part works is what lets you extend them and put statistical structure on top. Still, three questions remain open.

The first is too many variables. It is absolutely fundamental to linear regression that the sample size n exceeds the number of predictors p. This is not a case of breaking one small thing. When p > n, the matrix XᵀX cannot be inverted and the whole structure falls apart.

The second is predictive accuracy. Even when p < n, given an enormous set of features, is there a good subset of them to predict with? Can we improve the mean squared error by reducing variance?

The third is interpretability. Perfectly collinear features bring back the non-invertible matrix, but high collinearity is a problem too. The variance of the least-squares estimator is essentially

and with nearly collinear features, XᵀX becomes very small, so the variance explodes. Move the dataset a little and the estimate moves a lot.

These look like three separate questions, but the lecture argues they are all the same phenomenon: many features relative to the sample size. Sometimes that literally means p > n; sometimes it means that some of the features are very similar to one another. The area is also more modern than linear regression itself, because statistics and machine learning follow the data people want to analyse. From the 1990s, genetics brought datasets with perhaps 50 people in a study and thousands of gene-expression features for each one, and that renewed the interest in these questions beyond pure mathematics.

Following the textbook, the lecture looks at three families of methods:

1. Subset selection: identify a subset of the p predictors related to the response, and train only on those.

2. Shrinkage (or regularisation): train with all p predictors, but “shrink” the coefficients towards zero, which reduces variance and discards variables.

3. Dimension reduction: project the p predictors into a space of dimension M < p, and use those M projections as the variables of the model.

The first is a computational solution, with no mathematical guarantee and no mathematical sophistication either. The second is much more mathematical. It starts from the fact that least squares has no bias, since 𝔼[β̂] = β, but can have very large variance. Why not accept a little bias in exchange for a drastic reduction in variance? The third combines features rather than choosing among them. The book treats two methods of this kind, and this lecture covers the first.

Best subset selection: compare everything, and pay for it

We start with the computational family. The idea is silly but reasonable: compare all the models with k variables, for k = 0, 1, …, p, and keep the best. Once you fix k, the number of models to compare is

and you pick the one with the lowest training error, meaning the lowest residual sum of squares (RSS) or the highest R².

We usually distrust training error, but the reason is specific: it makes models with different numbers of features hard to compare, since the lowest training error always favours more features. Here the number of features is fixed. The risk of overfitting comes from comparing a five-feature model with a six-feature one and concluding that six is better. Among five-feature models, whichever does best even on training error probably has the most relevant features.

That settles the choice within each size, and leaves you with one best model per size. To choose among those, the argument no longer holds, because now the comparison is across sizes. You want to minimise validation error, not training error, otherwise you would always pick k = p.

K-fold cross-validation estimates this well, but you would have to run it for every one of those models, which can be heavy. Historically, though no longer today, this motivated alternatives for linear regression that estimate validation error directly. The first is Mallows’ Cₚ, or the Akaike Information Criterion (AIC):

The second is the Bayesian Information Criterion (BIC):

And the third is adjusted R²:

These come from asymptotic arguments with strong assumptions, such as Normal errors, and the derivations are not trivial. Deriving or even using them is not something that happens much. It is still useful to know what they are doing when someone mentions them: each one takes a notion of training error and reweights it by the number of features to make the comparison fairer. From this lecture’s point of view, they are simply ways of avoiding cross-validation when you have to.

The whole algorithm fits in three steps:

1. Let M₀ be the null model, with no predictors.

2. For k = 1, 2, …, p: train all the models with exactly k predictors and keep the one with the lowest training error, calling it Mₖ.

3. Among M₀, M₁, …, Mₚ, choose the model with the lowest validation error, using k-fold cross-validation, Cₚ (AIC), BIC or adjusted R².

There are two serious problems with this. The first is computational: searching every model means 2ᵖ possibilities, so unless you have something like p = 4, you are in trouble, and p = 100 or p = 1000 is out of the question.

The second problem is worse, because it is statistical. You have one fixed validation set, a random sample. If you test five, seven or ten models, the chance that one of them looks perfect on this sample purely by luck is very low. If you test thousands, that chance starts to grow: you try so many things that you will stumble on a model that happens to be very good on this one validation set and bad on everything else.

Testing too many models overfits the validation set itself, so a bigger search is not only slower but statistically riskier.

Forward and backward stepwise: searching a smaller space

Both problems suggest the same remedy: search fewer models, in a more directed way. Fewer models are cheaper to run and less likely to overfit. The tool is a greedy algorithm that adds or removes one feature at a time, and adding one at a time gives forward stepwise selection:

1. Let M₀ be the null model, with no predictors.

2. For k = 1, 2, …, p: consider all the models that extend the current one with one additional predictor, and keep the one with the lowest RSS (or the highest R²).

3. Among M₀, M₁, …, Mₚ, choose the model with the lowest validation error.

The change is in step 2. Once you have decided on your ideal model with three features, finding the one with four only requires looking at the features that remain. You do not reconsider every four-feature model. An example from the book shows the consequence:

Best subset can pick more than one new predictor at a time, while forward stepwise can only add one to what it already has. What is probably happening is that rating is highly correlated with the other features: it is very useful when you may only choose one, but once you may choose four, a better combination exists without it.

On training error you can guarantee something. Everything available to forward stepwise is available to best subset, so best subset never loses there. On test error you can guarantee almost nothing. Sometimes forward wins, sometimes best subset wins, and sometimes a sillier model beats both. On the test set, you run it and see what happens.

If you can walk forwards, you can walk backwards, and backward stepwise selection starts from the full model Mₚ. For k = p, p-1, …, 1, consider each of the k models that drop one predictor from the current one, and keep the one with the lowest RSS (or the highest R²). Then choose among the resulting models by validation error, as before.

Both directions choose among only

models, where the 1 is the null model. That is far smaller than 2ᵖ, so the risk of overfitting the validation set drops too, because the procedure is less flexible about which features it can combine. Neither guarantees the best model, though that may not be a problem. The two can also end up in different places: run both looking for the best five-feature model and you may get two different answers, because in a greedy algorithm the path matters.

There is one practical asymmetry. Backward selection cannot be used with p > n, because it begins by fitting the full model, which you cannot do. Forward selection can start small and stop at a model of size n, so for large p it is the only option.

Variations are easy to invent. Mixed stepwise may discard a selected predictor at each forward step, taking two steps forward and one back. Forward stagewise decorrelates the remaining predictors after adding one, before choosing the next. If some features are only useful together, you can add them as groups. A student suggested computing the correlation between features and keeping only one of each highly correlated pair, which the lecturer welcomed as domain knowledge, like blocks of genes that probably express the same thing. There is a catch, though: cross-validation assumes the validation data have never been seen, so computing those correlations on the whole dataset is a risk.

In this computational world you have few restrictions, but also few guarantees. Notice also what has and has not changed. The data, the class of linear functions and the squared-error loss are the same as before. The only change is how the minimisation is solved: sequentially, instead of in closed form, because the closed form needs an invertible XᵀX.

Ridge regression: change the loss, not the search

The stepwise methods still train many models, and they deal with the cost by shrinking the search. A recurring move in this course offers another way out, which is to change the loss function. The idea is to train one single model with all p predictors, with a loss that encourages many of the coefficients to be as close to zero as possible, and ideally exactly zero.

Why would that be a good idea? If the coefficient of a feature is zero, it is as if the feature did not exist. Tomorrow, to make a prediction, you would not even need to collect it. Forcing coefficients towards zero certainly introduces bias, because the expected value of the estimate is no longer β. But it can drastically reduce the variance of the estimators, and potentially lower the test error.

Ridge regression determines the coefficients by solving

The first term is the RSS of the linear model, exactly as before. The second term is new: λ times the squared L2 norm of the coefficients, which imposes a penalty in favour of shrinkage. Two extreme cases give the intuition. With λ = 0 the penalty disappears and you are back to ordinary least squares. With λ = ∞, any non-zero coefficient would make the loss infinite, so the only option is to set them all to zero, giving a constant prediction that ignores the features entirely. In between, λ is a hyperparameter that measures how much you care about fit versus shrinkage.

Remarkably, this problem has a closed form:

Before, the trouble was that XᵀX was not invertible when p > n. Simply adding λ along the diagonal, with λ > 0, makes the matrix positive definite and therefore invertible. After a long stretch of the course with no closed forms, in logistic regression and GLMs, we have one again, even when p > n. (Negative values of λ are a different matter: the lecturer called that case heretical and left it for after class.)

With λ close to zero, this is almost ordinary regression: limit sits above 400 and income near -300. As λ grows, the coefficients converge to zero, and past a certain point the prediction is essentially constant. The path need not be monotone, as rating shows. This is uncommon, and it echoes the earlier table: as limit is forced down, the information it carried is lost, and a previously less relevant variable can gain weight as the variables reconfigure. The only thing being imposed is that the L2 norm of the whole vector falls as λ grows.

How do we pick λ? The same way we pick any hyperparameter: try many values and use cross-validation.

Since the grid could have included values close to the extremes, an intermediate minimum suggests a model better than ordinary regression. This is where the opening questions start to combine: better prediction, a smaller effective set of features, and something easier to interpret.

As λ increases, bias rises and variance falls, so you are walking along the bias-variance trade-off. The lecturer added that you can prove there is always a non-negative λ giving a model strictly better on average than linear regression. The catch is that finding it is not trivial: you must estimate it by cross-validation, and the act of estimating it may be enough to lose the advantage. In that case, plain linear regression is worth keeping.

Two practical details

Ordinary least squares is invariant to scale. If you replace Xⱼ with cXⱼ, the optimal coefficient simply becomes β̂ⱼ/c, so metres or centimetres make no difference to the predictions. Ridge does not have this luxury. With a finite budget for the norm, you cannot multiply one coefficient by 100 and get away with it, so the same problem in metres and in centimetres gives different models. The fix is to put every predictor on the same scale, with empirical variance 1, before fitting:

This stops some coefficients being unfairly penalised more than others.

You may also have noticed that the penalty leaves out β₀. Nothing is mathematically wrong with penalising it, but the intercept acts as an estimate of the mean: if you had to predict a single number, it would be this one. In a genetics problem some features are relevant, but the intercept is always relevant, so it is not shrunk. The practical way to do this uses its least-squares form,

Substituting this into the residual of an observation gives

So you centre y and the columns of X, marked below with the subscript c, and solve

Then you add ȳ back when predicting. Any Python package already does this for you.

The lasso: coefficients that are exactly zero

Ridge solves a lot, but it never quite does what the shrinkage idea promised at the start: discarding features. The second shrinkage method changes one small detail. The lasso solves

This is almost the same problem, with the L1 norm in place of the squared L2 norm. Why bother with a separate name? Because from an applied point of view it is very different. Ridge shrinks all the coefficients jointly and smoothly, and essentially none of them is ever exactly zero. The lasso sets several of them to zero: not nearly zero, but exactly 0.0.

Given how this lecture began, that matters a lot. Suppose you run the lasso on a genetics dataset and find that many genes are irrelevant to, say, predicting cancer. Tomorrow you do not need to collect those features at all. With ridge that never happens, since every feature still matters, some of them only very slightly, and you cannot throw them out. The price is that the lasso, unlike ridge, has no closed form.

A student asked why these curves lose their smoothness. The answer is that the model is being forced to make a decision: include the variable or not. When a coefficient near zero becomes zero, that discrete decision has a systemic effect, and the model reconfigures in a non-smooth way. Ridge adjusts each coefficient smoothly, so nothing suddenly drops out. This ability to set coefficients to exactly zero is the main reason the lasso is so widely used.

L0, L1, L2: one problem written three ways

So far best subset and the penalised methods look like different worlds, one computational and the other mathematical. Rewriting the penalties makes them meet. Thinking in terms of Lagrange multipliers, which is where the λ comes from, ridge is equivalent to minimising the RSS subject to a budget s:

Here s is a hyperparameter of the same character as λ: for every λ there is an s giving an equivalent solution, and vice versa. The lasso is the same with the L1 norm:

And best subset selection, the purely computational problem we started with, turns out to fit the same form:

This last constraint is the so-called L0 norm. It does not penalise the magnitude of the vector but the inclusion of a feature: include it and pay one, leave it out and pay zero. So we are moving from L0 to L1 to L2. L2 gives a closed form. L1 does not, but it has an enormous advantage over best subset: the problem is convex, which means we know how to solve it efficiently.

One way to think of the lasso is as the convexification of best subset selection: a problem we would love to solve, replaced by a close one that we can.

The analogy with classification is direct. There, the 0–1 loss is what we want but cannot minimise, so we use the logistic loss. Here the L0 norm is what interests us, but we cannot optimise it, so we use L1. Computation, algorithms, optimisation and mathematics, which looked like totally different ideas, start talking to one another. Asked whether these methods recover the “right” coefficients of some generating process, the lecturer was careful: there is no such guarantee. What you get is a coefficient vector with a small norm that predicts very well.

Why does L1 produce zeros in the first place? The book justifies it with a picture.

Solving the problem means finding the best level curve that still touches the admissible region. Ridge’s ball is smooth, so short of very special circumstances, the curve touches it at a point where no coefficient is exactly zero. The lasso’s ball is a diamond, and its geometry pushes the best point to the extremes, which are exactly the points where a coefficient vanishes.

The lecturer treats this picture with some suspicion. The lasso is meant for very large p, thousands or millions of features, and it is not clear that the drawing looks like this in such dimensions. There are good mathematical reasons for the behaviour, and it is very clearly understood, but the visual intuition is a bit tricky. The picture does its job of showing that the two norms do very different things.

You can also place the norms on a continuum. L0 would be just two sticks along the axes, since a variable is either in or out. L1 inflates those sticks into the diamond. From there the ball gets rounder and rounder, and at the infinity norm it becomes a square. Asked why not use something between L0 and L1, the lecturer turned the question around: what algorithm would you use to optimise it? L1 is the trade-off, because it is convex.

Stepping back, the problem now is to minimise a penalised squared error,

with the L1 norm or, for the closed form, the squared L2 norm. It is convex, so it admits global solutions. Stepwise methods approximate best subset algorithmically; these approximate it by changing the loss. And nothing here is special to linear regression. A student asked whether GLMs have an equivalent, and they do: somewhere inside there is a linear form, so you penalise β there, and a convex function composed with a linear one stays convex. Predicting bicycle use across Rio de Janeiro with every feature you can collect, you would use L1 if you think few features matter, and L2 if you just want stability for better prediction.

Ridge or lasso? What each one assumes

If the two penalties are so close, how do you choose? Since the lasso forces some coefficients to be exactly zero, you would expect ridge to do better when most coefficients are non-zero. A simulation from the book confirms it.

In both methods, larger λ means more bias and less variance. Here the bias curves are almost the same, but ridge gains in variance by shrinking all the coefficients homogeneously. Now run the same simulation with only two non-zero coefficients, without knowing which two.

The picture reverses: both the bias and the variance of the lasso are lower, and the gain in error is large. The moral is that the mathematics is very similar but the assumption behind it is not. With ridge you assume many features are relevant, but shrinking them gives a better prediction. With the lasso you assume that only some are relevant and the others should be exactly zero.

Tuning works just as before, and in this example the model chosen by 10-fold cross-validation uses only two features.

A case you can solve by hand

Because the lasso usually lives in very high dimension, refined understanding of its mathematics is hard. The case n = p with X = I simplifies everything: multiplying by the identity just picks out each entry, so each coefficient can be minimised separately. Ridge then minimises

and the solution follows at once from the closed form, since XᵀX = I and Xᵀy = y:

The interesting part is that the lasso can also be solved here. It minimises

The absolute value has a kink at zero, with no derivative there. The technique is the subgradient, which amounts to pretending you can differentiate at zero, since any slope in the kink will do. Elsewhere the derivative is the sign function, and setting the first-order condition to zero gives three cases. For large positive values,

For large negative values,

And everything small collapses:

Here the two behaviours are plain to see. Ridge shrinks proportionally: the slope goes from 1 to 1/(1+λ), so a bigger yⱼ is pulled in by more. The lasso shrinks equally: anything very small collapses to exactly zero, and everything else moves by the same amount, whatever the size of yⱼ. This has everything to do with the L0 loss it approximates. In best subset, the size of a coefficient never mattered, only whether its feature was in. The lasso behaves similarly: some coefficients become zero, and the rest are not penalised more for being large.

Principal components regression: compress the features first

The third family takes a different route, somewhere between the other two. Instead of touching the search or the loss, it transforms X₁, …, Xₚ to improve the regression. Picture collecting everything you can about bicycle use in Rio: temperature, whether it is sunny, and much more. All those features look at the same phenomenon, which is how the weather is today. In genetics, many genes may relate to the same protein. Rather than choosing among them, why not combine them?

Suppose we build M < p new variables, each a linear combination of the original predictors, with constants φ still to be chosen:

Then fit an ordinary linear regression on these new variables, with coefficients θ:

If the constants are well chosen, this can beat the original regression. To see what is going on, substitute the definition of the z and group the terms by x:

So this is a special case of ordinary linear regression, where each coefficient is restricted to the form

The smaller M is, the tighter the restriction, which is natural when you compress everything into a couple of features. As with stepwise methods we are restricting a space, but it is no longer a computational one: it is the space of admissible parameter values, as if we were changing the class of functions the model can take. If the Zₘ summarise the predictors well, principal components regression (PCR) can reduce variance without costing much bias. Some bias does come in, as with ridge and the lasso, because β can no longer be anything. A student asked whether this also narrows confidence intervals. If the constants are chosen from the data, the lecturer answered, you cannot even talk about confidence intervals, because estimating in two stages loses the assumptions behind them.

The remaining question is how to choose these new variables. For now, the proposal is to do it without using Y, a decision to be revisited in the next lecture. One reason is to avoid overfitting, since Y is already used to estimate the θ. The choice is this: Z₁ is the direction that maximises the empirical variance of the data. The second direction also maximises the variance, subject to being orthogonal to the first, and so on.

Nearly all the variation in these data happens along the green line, while the blue direction looks like noise. Notice that nothing here says what we want to predict; we are only trying to summarise. The two variables carry practically the same information, so we could consolidate them into one.

Using a direction means projecting onto it: each point in two dimensions lands on the green line, which summarises it with one number. With only two dimensions, there is just one orthogonal direction left, but in higher dimensions you repeat the search many times and keep only the first M directions. The second component looks unimportant here, and plotting the original variables against each component confirms it.

Finding the directions, and choosing how many

The picture needs mathematics behind it. Centre the columns of X to mean zero, again marking the result with the subscript c. The empirical variance is then the matrix

The direction we want maximises the variance, among unit vectors:

The constraint is necessary, because otherwise scaling v would make the objective arbitrarily large. The solution is the first eigenvector of S, the one associated with the largest eigenvalue. If you have only taken one linear algebra course, you may think of eigenvectors through the characteristic polynomial, but here is another lovely characterisation: write the Lagrangian of this problem and what comes out is exactly the eigenvalue equation Av = λ v. Eigenvectors can be defined as the solutions of this optimisation problem.

The combination of the features with the largest variance is then

This really is a projection, much as linear regression is a projection on the column space, with the usual inverse term disappearing because v₁ has norm one. Requiring the next vector to be perpendicular to the first gives the second eigenvector, and so on. In matrix form, if the singular value decomposition is

then V holds the eigenvectors, and all the components come at once:

So the abstract-looking search for directions is nothing more than finding eigenvectors of the empirical covariance matrix, which is not hard at all with an SVD in hand. There is one last detail. As with ridge, scale matters: multiply ad spending by a trillion and both your eyes and the mathematics would suggest a different direction. So the columns are also normalised to variance 1, and the full algorithm reads:

1. Normalise the columns of X to mean zero and variance one.

2. Compute the SVD to obtain the eigenvectors in V.

3. For m = 1, …, M, obtain the m-th component by multiplying the normalised matrix by vₘ, the m-th column of V.

4. Return the predictions of the linear regression of y on the M components, with intercept ȳ:

How many components should you keep? Like λ, M is a hyperparameter, and cross-validation is the natural way to choose it. That is why the chapter on validation had to come before this one: the models now need hyperparameters, and they work only if you know how to choose them.

As you add components, the variance rises and the bias falls. And when M = p, you are essentially doing linear regression, so PCR is, once again, a way to trade bias against variance.

On the Credit data, the model with 10 components is the best, and 11 would be the usual regression. With 11 features, you can apparently reduce to 10 without losing any predictive performance.

The last slide compares the methods on a simulation built with n = 50 and p = 45, where y depends only on the first five principal components.

With PCR, adding components helps until you reach five. From there the bias stays flat, since y can already be expressed perfectly and extra features add nothing, while the variance keeps growing because the predictions depend more and more on the data. Ridge and the lasso are not the same figure, but they are very similar, and they tell the same story.

Conclusion

The lecture started from three questions that turned out to be one: what to do with many features relative to the sample size. Best subset answers it by brute force, comparing every model, at a cost that is both computational and statistical, and forward and backward stepwise make that search affordable by being greedy. Ridge and the lasso train a single model with a penalty: ridge has a closed form that works even with p > n, while the lasso sets coefficients to exactly zero and is the convexification of best subset. Which one wins depends on whether you believe many features matter or only a few. PCR compresses the features into directions of maximum variance, found as eigenvectors, and runs a regression on the first few.

For all their differences, these methods share a logic. Sometimes you restrict the optimiser, sometimes the loss, sometimes the class of functions, and each time you do it carefully, with a hyperparameter balancing the two sides, to trade bias against variance and beat plain linear regression. The lecturer finds this dance between the three objects very beautiful, not least because the ideas arrived at very different times, the lasso being relatively recent. What comes next is taking the linear model to its limit, followed by extensions in directions that are surprising and very good at prediction.

Thanks for reading! If you found this useful, consider following for more articles on machine learning and statistical modelling.

The GLM recipe goes beyond classification. Then: how to estimate test error with cross-validation, and uncertainty with the bootstrap.

Every image and example was taken from the following sources:

• Lecture 5 — GLMs e Métodos de Reamostragem — Machine Learning 2026 Summer Course at IMPA
• Machine Learning 2026, Prof. Paulo Orenstein

Introduction

These notes cover Lecture 5 of the IMPA Summer Course — Machine Learning 2026 (Paulo Orenstein), one article per lecture.

If you want extra context, the earlier articles in the series follow the course order:

1. Foundations of Statistical Learning: From Rules to Data, and What Makes an Estimate Good

2. Linear Regression: Least Squares, Inference, and How Far a “Simple” Model Stretches

3. When Linear Regression Breaks: Five Ways the Assumptions Fail, and How to Compute β̂ Without Inverting a Matrix

4. Classification Is Not More of the Same: Logistic Regression, Discriminant Analysis, and How to Judge a Classifier

Everything below is self-contained: you do not need to read the other posts to follow this one, as long as you are comfortable with likelihoods, the idea of training and test error, and a little probability.

The lecture has two halves. The first closes Chapter 4 of the course textbook with generalised linear models: the recipe used to build logistic regression, applied to count data, gives Poisson regression. The second covers all of Chapter 5, on resampling methods. It asks how to estimate test error and choose hyperparameters without cheating, moving through the validation set, leave-one-out and k-fold cross-validation, the one-standard-error rule, the wrong and right way to cross-validate, and the practical traps around all of it. It then asks how to measure an estimator’s uncertainty when no formula is available, which is the job of the bootstrap.

One recipe, many distributions: GLMs and Poisson regression

So far the course has sorted its methods by the kind of output. If y is a real number, we use regression; if it is 0/1 or one of a few categories, we use classification. What if it is neither? A very common case is a count: the number of films someone watches in a month, or the rounded number of hours a day they spend on some activity. Linear regression can predict negative or fractional counts, and there is essentially no way to stop it. Logistic regression is built for categories with no order, which a count is not.

Rather than invent a new method for every new type of data, the lecture proposes to generalise the recipe behind logistic regression. That recipe had four steps. Choose a distribution for Y given X (there, the Bernoulli). Write its log-likelihood. Make the term that multiplies y linear in the features. Minimise the negative log-likelihood. For counts, a natural distribution is the Poisson, which only takes the values 0, 1, 2, and so on:

Its parameter λ(x) is at once the mean and the variance of Y given X = x. That is new: in linear regression the mean and the variance were separate things, while here choosing the mean already fixes the variance. Taking logs gives the log-likelihood of one observation:

The term multiplying y is log λ(x), so that is what we make linear:

Equivalently, λ is the exponential of a linear predictor, which keeps it positive. Estimating λ now means estimating the coefficients. Maximum likelihood (with the log-factorial term dropped, since it does not involve the parameters) gives the objective

Once more, yᵢ multiplies something linear in the parameters, and once more the problem is convex in the coefficients, so gradient descent or Newton’s method converge. The lecturer called it a mathematical miracle the course would not explain: for almost any reasonable distribution, precisely those in an exponential family, the same recipe leads to a convex problem. If you do not know what an exponential family is, you can set up the problem and check convexity directly. With β̂ in hand you get λ̂, and from it any probability you like, such as the chance of zero events or of more than five. To estimate the probability that someone spends more than five hours a day on their phone, from features like how often they exercise or how many books they read, sum the fitted Poisson probabilities beyond five.

Compared with linear regression, three things change:

- Variance. Poisson regression is rigid: the conditional mean and variance are the same number, λ. Heteroscedasticity, which was a problem to defend against in linear regression, is here part of the design.

- Non-negativity. Negative values are never predicted.

- Interpretation. Increasing a feature by one multiplies the expected count by e raised to its coefficient, instead of adding the coefficient. That is simpler again than in logistic regression.

Bikes in the city

The textbook example is a bike-sharing system, like the public bikes on some streets in Rio. Whoever runs it needs to know how many users to expect at a station each hour, so bikes can be moved around. That means predicting integers: nobody will ever see 3.5 users.

Figure 1 shows why Poisson fits. At 2 a.m. there might be two users and almost no variation. At a busy hour, sun or rain moves the count a lot. The variability grows with the mean. The book fits a Poisson regression with dummy variables for each month and each hour of the day.

Negative coefficients are no problem; each one acts through e raised to its value. The data are from the US, so the warmer months are the busier ones, and the bikes are used far more during the day than before dawn (Figures 2 and 3). With a model like this, the operator does not have to wait for a dashboard to show that a station has run dry. They can predict how many users a station will see in two hours and move bikes in advance. They can also test hypotheses, for instance whether mornings are busier than afternoons.

What Chapter 4 was really about

The lecturer framed this as the twist at the end of the classification chapter. The methodology behind logistic regression was never a classification method. Pick a distribution for Y, write its log-likelihood, treat it as a loss and estimate the coefficients: that is a general philosophy. Interested in 0/1 outcomes, you chose the Bernoulli and got logistic regression. Interested in counts, you choose the Poisson and get Poisson regression.

Generalised linear models also come with confidence intervals and hypothesis tests, which are asymptotic, but work very well. There is no closed form for β̂. What you are guaranteed is that gradient descent, with a little care over the step size, reaches the right answer. Leaving linear regression can feel like losing things: first the formula, then the exact distribution of the estimator. But the world of models is enormous, and landing on something convex is a stroke of luck. The 0–1 loss was not convex, and the whole construction existed to swap it for a loss we can optimise. It is the course’s recurring theme again: the class of functions, the loss and the optimiser shape one another.

Asked whether the Poisson assumption can be tested, the lecturer said there are tests, but all of them require setting part of the data aside. In predictive work, unless you hold a strong prior, you run the model and see whether it works. Choosing Poisson should not be only about having the right support. It is also a question of whether the distribution suits how the data arise, just as Gaussian noise in linear regression was justified as the sum of many small independent effects. Asked about adding higher-order terms, the lecturer pointed to the usual risk: more terms, more overfitting, and you decide when to stop by looking at test error. That question leads straight into the rest of the lecture.

Why training error proves nothing

Chapter 5 changes hats. Until now the course has discussed methods: linear and logistic regression, LDA, QDA, naive Bayes. Now the topic is methodology: how to tell whether a method worked. In the lecturer’s experience, this matters much more than picking the right model. A poor model will be caught and discarded by a sound evaluation, but if the evaluation itself is wrong, everything is lost. The chapter has two parts, validation and the bootstrap, and both lean on computation rather than mathematics. As models grow more complex, expecting closed-form estimators and exact tests becomes naive.

The question left hanging at the end of the previous lecture motivates the first part. With so many classifiers available, how do you choose? The test set is meant to estimate performance in the real world. Say a doctor will use your model on patients only if its error is below 10%. If you use the test set to pick between LDA, QDA and naive Bayes, you keep whichever did best on it, and the error you report is optimistic: you are cheating. Choosing on the training set is worse still, because it pushes you towards far more predictors than you need.

Picture an intern who arrives with a fantastic prediction function. You will deploy it only if the error in production is below 5%, and the intern reports 1% on the training data. That number means nothing: a low training error is very easy to get, and you do not know what the intern did behind the scenes. The safe procedure is the one the course has used so far. Split the data at random into training and test sets, give the intern only the training part, and evaluate their model on the part you kept. Whatever they did — a neural network, a regression, overfitting or underfitting — the test estimate is fair, because they had no way to tune it.

Training error usually underestimates test error dramatically; the gap between the two is what tells you something.

The training error is not useless, though. Ask the intern for it alongside the test predictions, because the relationship between the two diagnoses underfitting and overfitting.

On the left of Figure 4 the test curve is still falling and stays close to the training curve, so there is no overfitting yet: high bias, low variance. On the right, the training curve keeps falling while the test curve climbs, which is the signature of a model that is too flexible. With complex models such as boosting or neural networks, looking at these two curves is the first diagnostic. A very low training error next to a high test error says to reduce flexibility.

Three sets, not two: the validation set

The intern also has decisions to make: logistic regression or naive Bayes, LDA or QDA, which polynomial degree. They should do what you did, at their own scale. Split the data they received into a training part and a validation part, fit on the first and compare candidates on the second. In practice the data end up in three sets:

- training data, used to fit the model for each candidate value of a hyperparameter;

- validation data, where the hyperparameter with the lowest error is chosen;

- test data, used once the final model is fixed, as the measure of future performance.

The distinction is between parameters, which the algorithm learns (the coefficients of a regression), and hyperparameters, which it does not (LDA versus QDA, the degree of a polynomial). Getting this wrong is, in the lecturer’s words, one of the most common and most basic errors in the field. Note also that each set holds less data than you have. The usual practice is to retrain the chosen model on everything before handing it to the doctor, since there is no reason to waste data, which makes the estimates along the way somewhat pessimistic.

A student asked how models trained with different losses can be compared. The evaluator fixes their loss, say mean squared error or the 0–1 loss, and announces it. Each colleague can train however they like, logistic loss included, but everyone delivers predictions, and those are scored with the evaluator’s fixed loss. The training loss is just one more choice inside each method.

The validation-set approach is the procedure used so far: split the data in two, train on the first part, compute the error on the second, and pick the candidate with the lowest error.

The lecture’s example is a polynomial regression of fuel efficiency (mpg) on horsepower in the textbook’s Auto data.

Figure 7 exposes the problem. A different random 10% held out gives a different curve. Every split agrees that degree 2 beats degree 1, but the level of error you would call acceptable changes a lot, so a random choice ends up driving your view of the model. As a student pointed out, one unlucky split can send a cluster of unusual points into validation that the model never learned to predict. The approach has three problems:

- the estimates vary a lot with the split;

- only part of the data is used to train;

- the validation error tends to overestimate the test error, because the model is trained on a fraction of the data and will later be retrained on all of it, where more data will make it better.

Leave-one-out and k-fold cross-validation

Instead of fixing one split forever, cross-validation rotates the validation role through the data. In leave-one-out cross-validation (LOOCV), for each point i from 1 to n, you train the model on every point except i, predict point i, record the error, and finally average the n errors.

Write ŷᵢ⁽⁻ⁱ⁾ for the prediction of point i from the model that never saw it. In regression the LOOCV estimate is the average squared error,

and in classification it is the average 0–1 loss,

LOOCV fixes two flaws of the validation set. It is not random: it always returns the same number. And it uses almost all the data for each fit, so it does not overestimate the test error; with n = 10,000, each model trains on 9,999 points.

It has two drawbacks of its own. The first is subtle. Averaging is supposed to reduce variance, but the n errors being averaged are not independent. The model trained on points 2 to n is almost identical to the one trained on point 1 and points 3 to n, so their errors are highly correlated, and averaging correlated quantities buys little. The second is blunt: LOOCV retrains the model n times. Some deep-learning models cost millions of dollars to train once, so that is out of the question.

Asked whether the gap between LOOCV and other estimates says something about generalisation, the lecturer said not as far as they knew. These are meta-methods for evaluating methods, the theory behind them is thin, and that is an open research problem. So there is no single answer about which to use. There are computational limits and there is practice.

The practical answer is the middle ground: k-fold cross-validation. Divide the data into k subsets, typically k = 5 or 10. For each fold, train on the other k − 1 folds and compute the error on the one left out, then average the k errors. Has anyone proved that 5 is the right number? The lecturer joked that we are born with five fingers on each hand and ten in total, and that this is not much of an exaggeration.

Every point is used for training at some stage, which the validation set never managed, and the cost is 5 or 10 fits instead of n.

Compared with Figure 7, the 10-fold curves in Figure 10 barely move from one split to another. Side by side, the two methods compare like this:

k-fold trains on less data than is available, so like the validation set it overestimates the error, though less badly. In LOOCV the training sets are so alike that averaging n folds reduces variance less than you would hope, whereas k-fold training sets overlap but differ by a non-negligible amount. On every count k-fold sits between the two extremes, and it is what people use in practice.

Choosing hyperparameters, and the one-standard-error rule

Do cross-validation estimates pick the same hyperparameter the test error would? The textbook checks with three simulated regression problems.

The curves in Figure 11 are not identical, but they are not far apart either. Since the point of validation is to choose a hyperparameter, here the flexibility, what matters is where each curve has its minimum. CV would pick roughly 6, 4 and 10 in the three panels; the cross on the blue curve shows what the test data would have picked, which you are not allowed to use. Not peeking costs something, but not much, and 10-fold and LOOCV make practically the same choices.

Classification behaves the same way (Figure 12): CV may pick an order of three where the test error would pick two, and choosing the minimum of the CV curve is a natural rule. There is a catch, though. 10-fold CV is itself random, and its minimum might sit where it does purely by chance. The fix is a heuristic — much of this lecture is, the lecturer admitted — called the one-standard-error rule.

With 10 folds you have ten error estimates, so besides their mean you can compute their standard deviation, which measures how much the error varies across folds; those are the bars in Figure 13. Taking the minimum would mean choosing p = 10 or so, but every model from 9 to 15 (honestly, up to 20) has essentially the same CV error, and one of them will come out lowest by luck. The rule is to find the minimum, agree that everything within one standard error of it is indistinguishable, and among those choose the simplest model. Here that is p = 9.

The wrong way and the right way to cross-validate

Cross-validation is a simple procedure, and that is exactly why it is easy to get wrong. The lecturer’s example comes from papers, published in journals such as Nature, that made this mistake.

You want to predict whether each of n = 50 people has cancer, using a nearest-neighbour classifier (one neighbour) and p = 5,000 gene-expression measurements. The proposal: using all the data, find the 100 genes most correlated with the outcome, then estimate the classifier’s error on those 100 genes with 10-fold CV. Is that reasonable?

Exaggerate the setup to see the problem. Suppose the outcome is decided by a coin flip and has nothing to do with genetics. All 5,000 features are then irrelevant, and any classifier built on them should be wrong about 50% of the time. The procedure above reports a 10-fold CV error of 3%. Imagine publishing a cancer predictor with that accuracy. It has nothing to do with the model; it is a fundamental error in the evaluation.

The selected genes are the 100 most correlated with the entire outcome vector, test folds included. Out of thousands, they are the ones that by chance looked most related to y, so the relationship is spurious, but every fold has been handed features that are good for its own targets. Whatever the test fold is, it no longer resembles a new sample, because the features were chosen using exactly its labels. The decision about which features to use was made after seeing the outcomes you are trying to predict.

The right way moves the selection inside the loop. Split the data into 10 folds. For each fold, using only the other nine, select the 100 best genes and train the classifier on them, then compute the error on the held-out fold. Average the 10 errors. Done this way, the estimate goes back to about 50%.

This feels odd, because different folds will choose different features. A student asked how a single final model could come out of that. The answer is that cross-validation evaluates a method, not a fitted model. If your method is “regression with forward selection of 100 features”, each fold runs that method and produces an error, and the average is the CV estimate for the method. To deliver a model, you run the same method on all the data. Whether it ends up with the same features as any fold is irrelevant; what was validated is the procedure. Which procedure is good for choosing features is the topic of the next lecture.

Any methodological choice learned from the data — such as which features to use — has to happen inside the cross-validation.

When the assumptions behind evaluation fail

Even a correctly run cross-validation rests on assumptions, and errors in evaluation are the critical ones. The lecturer did not go deep, but listed the main ways real projects break them:

- Data leakage: information about the label sitting in the features. The wrong-way CV above is one case. In a famous X-ray model that detected broken bones with startling accuracy, the model had learned to read annotations on the films, not the images.

- Sample bias: training and test data come from different distributions. In election polling some groups are more reluctant to answer. In one famous US election, supporters of a controversial candidate avoided admitting it on the phone; the polls pointed to the other candidate, and the result was the opposite.

- Selection bias: the training data are not a random sample. Paediatric studies cannot randomly assign babies to risks; nobody will take pacifiers away from some babies and give them to others to see whether they cause harm.

- Non-stationarity: the distribution changes over time — seasons, relevant politicians. With global warming, predicting the future from past data alone underestimates the temperature.

How do you catch these when the data are too abstract to spot an annotation? There is no universal answer. The course may give the impression that most of the work is clever models, but in a real project about 80% of the time goes into looking at the data. In one project segmenting fetal MRI scans, the model worked beautifully until it suddenly failed on a few cases. Piercings near the scanned region distorted the images and were very rare, and triplets were another rare case. Analysing the data first and then studying where the model’s errors are largest almost always reveals the problem. The field of interpretability, which tries to understand how complex models make their predictions, helps here too.

A related question came up about very imbalanced outcomes, such as a rare cancer: careless folds could put every positive case in one fold, and the model would never see one. There are two remedies. One is to genuinely randomise the order and trust each fold to look like the whole. The other is to stratify, forcing every fold to contain the same share of positives.

The bootstrap: uncertainty without formulas

The second half of Chapter 5 turns from test error to uncertainty. Until now, finding the standard error of an estimator was done case by case. Take the sample variance:

To get its standard error, the course assumed the data were Normal. Then a rescaled version of the estimator follows a chi-squared distribution with n − 1 degrees of freedom, and knowing the distribution answers every question: its variance, or the probability that it exceeds some threshold. But all of that rested on Normality. What if the data are not Normal, or the estimator is too complicated to have a known distribution, or has no closed form at all? An idea from the second half of the last century was to set the mathematics aside and use the computer instead; machines, the lecturer noted, have improved a little since the 1950s.

Two assets

The motivating example is from finance. X and Y are the returns of two assets.

You invest a fraction α in X and the rest in Y, for a return of αX + (1 − α)Y. You want a good return but also a safe one, so you choose α to minimise its variance. The variance is quadratic in α, and minimising it gives

You do not know these variances and the covariance, so you “put a hat on everything” and plug in the usual estimates:

Suppose your sample gives an estimate of 0.6. Can you trust it? With a different sample, would it move to 0.59, or to 0.3? With real money at stake, that is the question that matters. If you knew the joint distribution of X and Y, you could simulate many samples, compute the estimate on each, and look at how much it varies, which is what the textbook can do because its data are simulated.

In real life you have one sample and no distribution.

Pulling yourself up by your bootstraps

The bootstrap’s answer is to resample from the original sample and pretend that this is equivalent to drawing new samples. You have n points, so you draw n points from them with replacement. Some points repeat and others are left out. You compute the estimate on this new set, repeat the whole thing m times, and use the standard deviation of the m estimates as the standard error. In effect you simulate the uncertainty without ever knowing the estimator’s distribution.

The name, as the lecturer told it, comes from the stories of Baron Munchausen. Stuck in quicksand and about to drown, the baron pulled himself out by his own bootstraps. The method seems to do something similar, creating a sense of variability out of a single sample.

Three observations are far too few in practice, but they fit on a slide (Figure 16). Each resampled table is treated as if it came from the same distribution as the original, each gives an estimate, and the spread of those estimates is the uncertainty.

Why the bootstrap works, and when it doesn’t

Does it work? In the textbook example, with m = 1,000 resamples (a usual number), the average of the bootstrap estimates is 0.5996, close to the true value of 0.6. The standard deviation of the estimates is 0.083, and that is the bootstrap standard error:

That is close to the true value. More useful still, the whole distribution of the bootstrap estimates approximates the true sampling distribution of the estimator.

The secret is not only that the standard error comes out right. In many cases the resampled data sets are distributed much like the real data, and when the distributions are close, any statistic you extract from them is close too (Figure 17).

Figure 18 puts this more formally. In the real world, an unknown population generates your data and you compute an estimate, and not knowing that population has been the problem all course long. In the bootstrap world, the role of the population is played by the empirical distribution of your sample. Treat the sample as if it were the whole population, like a histogram: if the value 3 appeared 20 times and the value 4 appeared 100 times, draw according to those frequencies. Sampling n points from it independently is what produces the repetitions. If the estimated population resembles the real one, statistics computed in the bootstrap world resemble the real ones.

Are there guarantees? There are, but they belong to a statistics course. Even stating what kind of convergence is meant is not trivial, since something empirical is converging. In practice, the lecturer said, it almost always works, with an asterisk discussed below. As for how many resamples, the limit is computing power: 1,000 is common and 10,000 is better if you can afford it. The more the better.

The payoff is large. All the work of deriving t and F tests, and the tests for logistic regression the course never derived, existed to find the distribution of a statistic, and the bootstrap now provides that distribution by computation. That will matter more and more as models become too complex for formulas. The bootstrap can give prediction intervals too, though more modern methods exist for that. It was introduced alongside a close relative, the jackknife, named after the pocket knife because it is useful for so many things.

The asterisk

The bootstrap is a non-parametric way to get the distributional properties of an estimator — standard errors, confidence intervals and the like. It relies on the resamples behaving like genuine independent samples, as if you could keep asking the doctor for one more data set. That approximation usually holds, but it can fail:

- Small samples. With n = 3, as on the slide, the resamples cannot capture the variability.

- Dependent data. Resampling points independently destroys the dependence. In a series where the temperature is rising, drawing points at random scrambles the temporal structure, which is the most important thing in it.

- Skewed, heavy-tailed distributions or rare events. The bootstrap distribution may not reflect them. The convergence results behind the method need some regularity; in practice this is rarely a problem unless the distribution is truly strange.

- Discontinuous statistics, such as the sample minimum or maximum. The maximum of a resample rarely matches the maximum of the full data, so the bootstrap underestimates it. For means, standard deviations or the probability of exceeding a threshold, it generally works.

There are fixes, including the jackknife, the parametric bootstrap and the block bootstrap for dependent data.

Instead of resampling points, the block bootstrap resamples whole blocks, hoping each block preserves part of the time dependence. In Figure 19 blocks 1 (green), 2 (blue), 3 (red), 4 (pink) and 5 (yellow) are cut and pasted whole. It is not perfect: the right block length depends on how long the dependence lasts, which nobody knows in advance, and smarter variants exist. For the independent data this course works with, the plain bootstrap is enough.

Conclusion

The lecture began by showing that logistic regression was never just a classification method. Choose a distribution for the response, write its log-likelihood, make the term multiplying y linear and minimise: the result is a generalised linear model. With the Poisson distribution, that recipe produces a convex, interpretable model for counts whose variance grows with its mean by design.

The rest of the lecture was about not fooling yourself. Training error proves nothing; a held-out test set is the only fair judge, and the gap between the two diagnoses under- and overfitting. Choosing models needs a third set, and a single validation split is noisy and pessimistic. Leave-one-out is deterministic but expensive and correlated, and k-fold cross-validation is the practical compromise, best combined with the one-standard-error rule. Every data-driven step, feature selection included, must happen inside the folds, otherwise a coin flip can look like a 3% error rate. Even then, leakage, biased samples and a changing world can undo an evaluation that looks correct. For uncertainty, the bootstrap replaces derivations with resampling: treat the sample as the population, recompute the estimator many times, and read the uncertainty off the spread, while remembering that it can fail for tiny samples, dependent data, extreme distributions and discontinuous statistics.

Both tools are computational and need few assumptions, which is why they will carry over to the far more complex models ahead. They are also how you protect yourself when a team builds models for you: a set nobody has seen is exactly how data competitions score their entrants. The lecturer recalled weather-forecasting contests where you predict a month ahead, impossible to game, and one where you had to “predict the past”, which was horrible for everyone because cheating was so easy. The next step is the question that has been hovering since feature selection came up: when there are many features, how do you choose which ones to use? With cross-validation in hand, the course is ready to answer it.

Thanks for reading! If you found this useful, consider following for more articles on machine learning and statistical modelling.

The loss classification wants can’t be optimised. This lecture rebuilds it around probabilities and asks how to judge a classifier.

Every image and example was taken from the following sources:

• Lecture 4— Métodos Lineares para Classificação— Machine Learning 2026 Summer Course at IMPA
• Machine Learning 2026, Prof. Paulo Orenstein

Introduction

These notes cover Lecture 4 of the IMPA Summer Course — Machine Learning 2026 (Paulo Orenstein), one article per lecture.

If you want extra context, the earlier articles in the series follow the course order:

1. Foundations of Statistical Learning: From Rules to Data, and What Makes an Estimate Good

2. Linear Regression: Least Squares, Inference, and How Far a “Simple” Model Stretches

3. When Linear Regression Breaks: Five Ways the Assumptions Fail, and How to Compute β̂ Without Inverting a Matrix

Everything below is self-contained: you do not need to read the other posts to follow this one, as long as you are comfortable with least squares, logarithms and exponentials, gradients, and a multivariate Normal density.

The first week of the course was about regression. This lecture, which follows Chapter 4 of the course textbook, turns to classification and argues that it is not “more of the same.” The natural loss cannot be optimised directly, and bolting linear regression onto 0/1 labels goes wrong. That leads to modelling the probability of each class instead, either directly, with logistic regression and its multi-class version (the discriminative methods), or through Bayes’ rule, with LDA, QDA and naive Bayes (the generative methods). We then compare them all and finish with how to evaluate a classifier when a single error rate hides most of the story. Generalised linear models, also on the lecture’s agenda, were deferred to the next class.

Why classification is not “more of the same”

In classification we predict a qualitative variable: a class that may have no natural order at all, like eye colour taking the values brown, blue or green. Code them as 0, 1 and 2, and the job looks familiar: find a function f̂ that takes features X and returns one of those three codes.

Returning a bare label throws information away, though: it says nothing about how confident the method is. “Your eyes are blue, and I am 90% sure” and “your eyes are blue, and I am 55% sure” are very different statements. So we usually predict a probability of belonging to each class — a real number between 0 and 1, smoother and more informative than a label. There is a second reason, which turns out to matter more: a discrete function is miserable to optimise, and probabilities smooth the problem out.

Figure 1 shows the idea on the textbook’s credit-card data.

The setup keeps the three decisions of every method in the course — a class of functions, a loss and an optimiser — but the predictive function gains a two-level structure:

What we really model is p̂(x), the estimated probability of class 1, and the prediction is the most probable class. We would like p̂ to be as close to linear as possible, but it cannot be exactly linear: probabilities cannot be negative or exceed 1, and they must sum to 1. A function h, still to be chosen, will enforce that.

The natural loss is the 0–1 loss: pay 0 for the right class and 1 for the wrong one. The obvious plan is to find the coefficients that minimise its average over the sample, the way we did with squared error. It does not work. Figure 2 shows why.

The loss is neither convex nor smooth, and its derivative with respect to p̂(x) is zero almost everywhere. By the chain rule, the gradient with respect to any coefficient is that zero times something, so it vanishes almost everywhere no matter how p̂ is parametrised. Gradient descent has nowhere to go, and the comfortable world of linear regression — convex, unique solution, closed form — is gone.

The way out is to rethink the loss, using the interplay between function class, optimiser and loss that the course keeps returning to. We train on a surrogate loss: something close to what we want but mathematically better behaved.

The obvious surrogate, and why it fails

The first candidate is the only loss we know so far, the squared loss. Run a linear regression with Y ∈ {0, 1} and use the fitted values as p̂. It is even aligned with the goal, because for a 0/1 variable the conditional expectation is exactly the probability of class 1:

Then threshold at 0.5. The lecture lists three problems with this. There is no obvious way to extend it to more categories, such as {stroke, cancer, diabetes}. The fitted values need not lie in [0, 1], so they are not probabilities: one negative coefficient and a large positive feature can produce a “probability” of −8. And the loss punishes the wrong things. If the true class is 1 and the model says 3 — so confident that it overshot the ceiling — it pays (3 − 1)² = 4. Had it predicted 0, which is much worse in any sensible sense, it would have paid only 1. Confident correct answers get over-penalised.

The underlying issue is that we are forcing a method built for real-valued y onto a problem that is intrinsically different. We need a different loss and, with it, a different predictive function.

From a coin flip to the logistic loss

The lecture’s proposal is to think statistically about what Y is. A variable that only takes the values 0 and 1 has essentially one possible distribution, the Bernoulli: a coin flip that lands on 1 with probability p(x) and on 0 with probability 1 − p(x). So we write

There is not much choice here: if Y is binary, it is Bernoulli. What we get to model is the function p(x).

This suggests a loss. Maximum likelihood picks the p̂ that makes the observed sample as probable as possible. The sample is the only thing we know actually happened, so we make it as unsurprising as we can. The probability of observing a particular label y can be written compactly as

which is just a tidy way of saying “p(x) if y = 1, and 1 − p(x) if y = 0”, since one of the two exponents is always zero. We minimise throughout the course, so we take the negative. We also take the logarithm, which is monotone and so costs nothing. The result is the loss

the logistic loss. It is a smooth substitute for the 0–1 loss, and because it is smooth we can optimise it. It is also more informative. It distinguishes being wrong with conviction from being wrong while unsure, which the 0–1 loss, recording only hit or miss, cannot do. Figure 3 makes this concrete.

The small price paid when right but hesitant is, in the lecturer’s words, the problem telling you there is still something to learn. This is the first big break from linear regression, which kept its natural loss to the end; here a discrete loss is deliberately traded for a smooth one.

Making the log-odds linear

Averaged over the sample, and with the terms regrouped so that what multiplies yᵢ stands apart, the objective is

Minimising over every possible function is hopeless, so, as in linear regression, we impose linear structure. We do not make p itself linear, because that brings back negative probabilities. We make linear the term that multiplies y, the logarithm of the odds:

The lecturer was candid that the deep reasons for this choice involve exponential families and natural parameters, material for a full statistics course. One consequence can be checked by hand, though. Write s for the linear predictor, exponentiate both sides, and collect the p̂ terms:

Therefore

This is the logistic function applied to the linear predictor, and it is exactly what a probability should be: positive, below 1 because the denominator exceeds the numerator, and summing to 1 with the probability of class 0, which is one over the same denominator. If you would rather skip the statistical mysticism, read it as a normalisation that squeezes a linear score into [0, 1].

Substituting this form back into the loss turns a search over functions into a search over coefficients:

Here βᵀxᵢ is shorthand for β₀ + β₁ xᵢ₁ + ⋯ + βₚ xᵢₚ. And then a gift arrives: this problem is convex. We never asked for convexity; it was earned the moment we chose to model linearly the term that multiplies y, and it will reappear every time that recipe is followed. Modelling another parameter, such as the mean, would instead impose constraints like keeping p̂ in [0, 1] and could make the optimisation non-convex.

Accuracy is the loss we want but cannot optimise; the logistic loss is the smooth version we can.

Fitting it: convexity, gradient descent, and Newton in disguise

With a convex objective, we are back near linear-regression territory, minus the closed form. Here is the loss again:

Gradient descent only needs the first derivative. The second derivative is what proves convexity. If vector derivatives feel unfamiliar, differentiate with respect to each βⱼ separately and stack the results; it comes to the same thing.

The first term differentiates to yᵢ xᵢ. For the second, the chain rule gives one over one-plus-exponential, times the exponential, times xᵢ, and that ratio of exponentials is precisely p̂(xᵢ). So the gradient is a sum of residuals weighted by the features:

Here X is the feature matrix, y the vector of labels and p̂ the vector of fitted probabilities. Differentiating once more, only p̂(xᵢ) depends on β, and the derivative of the logistic function turns out to be p̂(1 − p̂). The Hessian is

where W is the diagonal matrix whose entries are p̂(xᵢ)(1 − p̂(xᵢ)). Those entries are positive, so the Hessian is positive semidefinite, and the loss is convex in β, as promised.

Gradient descent, a first-order method, starts anywhere and repeatedly steps against the gradient with a step size ηₜ:

Convexity guarantees it converges. Newton’s method, a second-order method, also uses curvature: it steps by the inverse Hessian times the gradient. Plugging in the two expressions above, the factors of 1/n cancel:

Now comes the part the lecturer clearly enjoyed. Write the current β̂ₜ as (XᵀWX)⁻¹XᵀWXβ̂ₜ, which changes nothing. Also slip W W⁻¹ in front of the residual, which is allowed because W is diagonal with positive entries. Both terms then share the factor (XᵀWX)⁻¹XᵀW:

If W were the identity, this would read (XᵀX)⁻¹Xᵀzₜ: a plain linear regression of zₜ on X. With a diagonal W it is a linear regression with weights, the solution of a weighted least-squares problem:

The i-th observation effectively enters with a weight equal to the square root of p̂ᵢ (1 − p̂ᵢ). This is iteratively reweighted least squares (IRLS). Buried inside a classification problem, every Newton step is a least-squares fit, which we know how to solve efficiently without inverting a matrix — with the QR or SVD decompositions, or with gradient descent. Fitting a logistic regression with Newton’s method is walking through parameter space solving one linear regression per step.

Newton’s method for logistic regression is a sequence of weighted linear regressions.

The price relative to linear regression is that the closed form (XᵀX)⁻¹Xᵀy is gone, and with it the exact distribution of β̂. All we have is an iterative method, run until it converges.

Reading the coefficients, and why “student” flips sign

Run logistic regression on the credit-card data in Python or R, modelling default as a Bernoulli whose probability is the logistic function of an intercept plus a coefficient times balance, and you get a familiar-looking table.

The standard errors and p-values come from asymptotic tests, not derived in class: they generalise the linear-regression calculations and again trace back to modelling linearly what multiplies y. The catch is the word asymptotic: the test is exact only as n goes to infinity, where the estimator’s distribution becomes Gaussian. It is an approximation, but one scientists use routinely, especially with large samples.

To read the table, compute probabilities. At a balance of 1,000 dollars, plugging the coefficients into the logistic function gives

and the same calculation at 2,000 dollars gives 0.586. Doubling the balance takes the probability of default from under 1% to over 50%.

Now fit a different model, with a single 0/1 feature: whether the person is a student.

A student gets a probability of default of 0.0431 and a non-student 0.0292 — “slightly sad news,” the lecturer told a room full of students. But these are two separate regressions. Using balance and student status together needs a multiple regression, fitted exactly as before by maximum likelihood.

With balance and income in the model, the student coefficient is negative. A bank that already knows your balance and income should be more willing to lend once it learns you are a student. The explanation is a confounding variable, shown in Figure 4.

Three facts combine here. Students carry larger balances: they are starting out, earn little and pay for their education, “in particular in the United States.” People with large balances default more. Yet among people with the same balance, students default less — of all the people who owe money, in the lecturer’s words, students are the best: young, starting to earn, investing in an education that will pay off.

Alone, being a student goes with owing money, and owing money goes with default, hence the positive coefficient. With balance in the model, the coefficient is read conditionally on the other features: at a fixed balance of, say, 1,000 dollars, you would rather lend to a student. The confusion comes from the correlation between the features.

The lecture closed this part with the method’s weak spots:

- Collinearity inflates the variance of the coefficients, as in linear regression; with perfect collinearity there is no solution.

- Perfect separation is a new problem. If the classes can be split perfectly, infinitely many curves reach zero training error, the slope wants to go to infinity, and gradient descent keeps stepping without ever converging. This is always the case when the number of features is at least the number of samples.

- Interpretation is harder than in linear regression. A unit increase in a feature no longer adds a fixed amount; its effect on p̂ passes through an exponential over one plus an exponential.

- Inflexibility: for some problems, such as classifying handwritten digits (MNIST, later in the course), the model hits a performance ceiling.

More than two classes: multinomial regression

What if there are more than two classes, say emergencies that are heart attacks, overdoses or epileptic seizures? This was one of the failures of the regression shortcut, and here the statistical framing pays off: we just need a distribution with support on K classes. That is the Multinomial with a single draw, in place of the Bernoulli, with one probability per class, all summing to 1:

Only the factor of the observed class survives. Once the first K − 1 probabilities are known, the last is one minus the rest, so we single out class K and model linearly the log of each class’s probability relative to it:

Each of the K − 1 classes gets its own p + 1 coefficients, for (p + 1)(K − 1) parameters in total. This is exactly what logistic regression did with two classes: there, class 0 was “one minus class 1”, and that is where the “1 −” in the log-odds was hiding. Solving for the probabilities in the same way gives

for k = 1, …, K − 1, while the last class gets 1 over the same denominator. Class k uses its own coefficients, and the denominator sums over all classes. The probabilities lie in [0, 1] and add up to 1. The negative log-likelihood has the same shape as before — a linear part plus the log of one plus a sum of exponentials — the parameters come from maximum likelihood with gradient descent or Newton, and the problem is still convex.

Where the regression shortcut had no answer for many classes, the generalisation here came for free: put a sensible distribution on Y and let it suggest the functional form. The lecturer called logistic regression perhaps one of the most important methods in machine learning, one the course will keep extending, all the way to neural networks.

A different route: generative models and LDA

Everything so far estimated ℙ[Y ∣ X] directly, by assuming a functional form for it. The second half of the lecture takes a route with a very different structure. Estimate two other things instead:

1. ℙ[X ∣ Y]: given the class, how are the features distributed?

2. ℙ[Y]: how likely is each class before seeing any features?

Then combine them with Bayes’ rule:

The first ingredient is the strange one. Until now the features were fixed and the randomness lived in Y. Purely for modelling purposes, we flip that and treat X as random given the class. Take predicting whether a patient develops diabetes next year from their insulin level and how much they exercise. The usual question is: given this insulin level and no exercise, what is the probability of diabetes? The flipped question is: given that this person will develop diabetes, what is their insulin level likely to be today, and do they exercise?

The denominator plays no role in the decision. We want the most probable class for a given x, and the denominator is the same for every class. Because these methods model how the features are distributed within each class — enough, in principle, to generate features — they are called generative. Methods that model ℙ[Y ∣ X] directly, like logistic regression, are discriminative.

Linear discriminant analysis

The first generative method, linear discriminant analysis (LDA), makes two assumptions. The features within each class follow a multivariate Normal distribution, and the probability of each class, πₖ, is estimated by its share of the training data. If 90% of a representative sample is not diabetic, then with no features at all the best guess for a random person is a 10% chance of diabetes.

A student asked whether anything was being assumed about the spread of the Gaussians. Something is, and it is the heart of the method: all classes share the same covariance matrix Σ. The ellipses in Figure 5 have the same shape and differ only in where they sit. The class therefore enters the density only through its mean μₖ:

By Bayes’ rule, the probability of class k given x is a constant times πₖ fₖ(x). Now absorb into the constant everything that does not depend on k: the power of 2π in front, the determinant of the shared Σ, and, after taking logs and expanding the quadratic form, the term in xᵀΣ⁻¹ x. What remains is the discriminant function:

Predict the class with the largest δₖ(x). Another student asked how terms that depend on x can be thrown away. The answer is that for a given x we only compare classes, so a term added equally to every class cannot change the winner. If two discriminants are 1.2 and 1.5, adding 10 to both still picks the second.

Where does “linear” come from? Look at the points where you are indifferent between classes k and l, that is, where δₖ(x) = δₗ(x). On both sides x appears only multiplied by a fixed vector, so the set of such points is a hyperplane. The indecision between green and yellow in Figure 5 is exactly a straight line. Gaussian classes with a shared covariance produce linear decision boundaries.

To use it we estimate three things from the sample:

These are the class’s share of the data and its centre of mass. The covariance is the unbiased estimator, pooled across classes because they share it, with each point centred on its own class mean:

Then compute every class’s discriminant at x with these estimates plugged in and pick the largest. Both the modelling and the reasoning are a world away from logistic regression.

When to prefer LDA

The lecture gives two situations. First, when the classes are well separated, LDA is more stable than logistic regression. Perfect separation is exactly what breaks logistic regression, whereas for LDA it is the easy case: put one Gaussian in each corner. Second, with a small sample and features that are roughly Normal within each class, LDA does better.

The second point reflects a philosophy that runs through the course: strong parametric assumptions beat assumption-free methods when they are reasonable and lose when they are wrong. LDA needs a strong one, because it models a p-dimensional X given Y — far more to estimate than logistic regression’s one-dimensional Y given X — and Gaussianity with a shared covariance keeps that manageable. Asked whether the assumptions could be tested afterwards, the lecturer gave a machine-learning answer: compare the model with others and judge by its predictions.

Loosening LDA: QDA and naive Bayes

LDA’s recipe is “class prior times class density, normalised,” with Normal densities sharing Σ. Changing the shape of the densities generalises it in two directions.

Quadratic discriminant analysis (QDA) gives each class its own covariance Σₖ. That is the only change, but it matters. The quadratic term in x and the log-determinant can no longer be absorbed into a constant, because both now depend on the class:

The boundaries become quadratic. Figure 6 shows what that buys and what it costs.

QDA gains flexibility at the price of many more parameters. A student asked why, in the left panel, QDA does not simply estimate near-identical covariances and recover the line. The lecturer compared it to polynomial regression: a quadratic fit could produce a straight line, but given the extra flexibility it will use it to chase the data, which is overfitting to some degree.

Naive Bayes takes a different shortcut. The real burden of discriminant analysis is a distribution over a p-dimensional vector: few are tractable (the Gaussian is one, hence its use), and you must estimate every correlation between features — insulin with glucose, exercise with age — a number of parameters quadratic in p, per class for QDA. Naive Bayes assumes the features are independent within each class, so the joint density factorises into one density per feature:

That is useful when p is large. With Normal marginals, independence means a diagonal covariance, and the discriminant becomes a sum of one term per feature:

up to a constant. But the Normal case is where you least need naive Bayes, since LDA and QDA exist. Its real strength is that each marginal can be anything. Fitting a Gaussian to a 0/1 variable like “exercises or not” makes no sense; naive Bayes can use a Bernoulli for it and Gaussians for the rest, or estimate a marginal with a histogram. That suits dummy variables and text, where predicting the next character is inherently discrete. The independence assumption is strong, but in practice it usually works — and, asked how to test it, the lecturer again answered: use a test set and see whether the predictions are good.

Comparing the methods: four ways to write the log-odds

We now have four classifiers built on quite different ideas. The lecture compares them with one rewrite. Choosing the class that maximises ℙ[Y = k ∣ X = x] is the same as maximising its log-ratio against a reference class K, so look at what each method assumes that log-odds looks like.

For logistic regression it is linear, by construction:

For LDA, the ratio is class k’s prior times its Gaussian, over the same product for class K, with the same Σ in both:

A ratio of exponentials is the exponential of the difference. The term xᵀΣ⁻¹ x appears in both the numerator and the denominator and cancels, leaving only a part that is linear in x. So LDA’s log-odds is also linear. Logistic regression and LDA have the same flexibility — what one can represent, the other can too — and differ only in how they fit the parameters to data. Nothing in the derivations so far made that obvious.

For QDA the same calculation leaves the quadratic terms standing, because the covariances of class k and class K differ:

For naive Bayes, the log of a product becomes a sum, one function per feature:

Each gₖⱼ depends on which marginals were chosen, so it can be linear, quadratic or something else. But the features act on the log-odds independently: there is no product of xⱼ with xⱼ₊₁, so the features never “talk to each other,” whereas QDA’s quadratic terms let them.

Put together, these four forms explain how the methods are related:

- Logistic regression is not generative. It never models X given Y.

- LDA is a special case of QDA (QDA could, improbably, estimate the same covariance for every class), but it can perform better. If the classes really do share a covariance, the lecturer would bet on LDA. QDA spends part of its sample estimating flexibility it does not need, which is the intuition for why a more flexible method can lose to a more rigid one.

- Every classifier with linear boundaries is a naive Bayes, LDA included: choose each gₖⱼ linear and the naive Bayes form reproduces the linear one. With Normal marginals that share a variance per feature across classes, naive Bayes is LDA with a diagonal covariance, since then the two make the same assumptions.

- QDA and naive Bayes are not special cases of each other. QDA has interactions; naive Bayes has arbitrary per-feature functions.

The methods make different assumptions, and none of them is uniformly better.

What the experiments say

The textbook makes this concrete with six simulated scenarios, shown in Figure 7.

The lecture did not go through the scenarios one by one. The point is that when you change how the data are generated so that they match one method’s assumptions, that method wins. There is no universally best classifier; there are assumptions, and a method is good when its assumptions are.

A student objected that with enough data the most flexible method should win. The answer was a qualified yes. An unnecessarily flexible method spends statistical information on parameters it did not need, so it loses to a simpler method with the right assumptions — but as the sample grows, that gap shrinks. That, the lecturer said, is one of the great revolutions of recent years: with datasets the size of the internet, flexible methods probably beat very parametric ones, which is why neural networks with billions of parameters sit behind modern AI, against the tens or hundreds in the textbook’s examples. In this lecture’s problems, the right parametric assumptions still pay off — a fundamentally empirical intuition, the lecturer added, since for neural networks there is very little theory.

The moral is a warning against downloading several models, keeping the best score and shipping it. Understand what each method assumes, look at your data, and match the model to it.

Judging a classifier: confusion, thresholds, and the ROC curve

Every method above was trained on a surrogate — the logistic loss, or a Gaussian likelihood — because the loss we care about could not be optimised. To evaluate, we go back to that loss, the 0–1 error rate: the fraction of observations whose predicted class differs from the true one.

It describes what we want, but it is not very informative when the errors are unbalanced between the classes. Take a real IMPA problem: detecting deforestation by comparing this month’s sub-images of Brazil’s forests with last month’s. Nearly every sub-image shows no clearing, so a model that always answers “no deforestation” is right ninety-something percent of the time — excellent on paper, and useless. Likewise, almost everybody repays their loans, so a bank model that lends to everyone scores nearly as well as one that took months to build, and the 0–1 loss cannot tell them apart.

A better description is the confusion matrix, which separates the two ways to be right from the two ways to be wrong.

The 0–1 loss weighs a false positive and a false negative the same, and in real problems we often care much more about one than the other.

An example with LDA

Use LDA to predict default for the 10,000 people in the textbook’s data, predicting “yes” when the estimated probability of default exceeds 0.5. There are 333 defaulters (positives) and 9,667 non-defaulters (negatives). The training error is low: 23 false positives plus 252 false negatives, out of 10,000, is 2.75%. But a model that always says “no” errs only on the 333 defaulters, about 3.33%, so a model that ignores the features is almost as good. The source of the error is the false-negative rate: of the 333 people who defaulted, the model missed 252, about 76%. It largely fails to recognise the case a bank cares about most.

The cheap fix is to lower the threshold — call someone a likely defaulter when the estimated probability exceeds 0.2 instead of 0.5. This happens without touching the model. The probabilities p̂ are already estimated; only the rule that turns them into a class changes.

The false-negative rate drops from about 76% to 138 out of 333, about 41%. The price is paid elsewhere: the false-positive rate rises from 23 out of 9,667 (about 0.2%) to 235 out of 9,667 (about 2.4%). The overall error rate, counting the off-diagonal cells, rises to 235 plus 138 out of 10,000, about 3.7%. That is expected, since the most-probable-class rule at 0.5 is what minimises the 0–1 loss. Forcing more positive predictions trades errors on one side for errors on the other. Whether that trade is worth making depends on the problem. A court, for example, operates under a presumption of innocence: our society decided that jailing an innocent person is much worse than freeing a guilty one, so the two errors are deliberately weighed differently.

The whole range of thresholds

Figure 8 sweeps the threshold from 0 to 0.5.

At a threshold of 0, essentially everyone has a probability above it and is predicted positive, so there are no false negatives at all. As the threshold rises, the model becomes less eager to say “positive”: false negatives climb and false positives fall. The overall 0–1 error combines the two.

The ROC curve packs every threshold into one picture. It plots the true-positive rate against the false-positive rate as the threshold moves (ROC stands for receiver operating characteristic, a name from engineering that, as the lecturer noted, has nothing to do with anything here).

At threshold 0 everything is predicted positive, so both rates are 1 and the curve starts at the top-right corner. At threshold 1 nothing is predicted positive, both rates are 0, and it ends at the bottom-left. In between, the closer the curve hugs the top-left corner, the better. That corner means catching every positive without a single false positive: perfect prediction.

Until now every evaluation needed a threshold. But if someone hands you a classifier and asks which threshold to use, there is no good answer without knowing how they weigh false positives against false negatives — and whoever builds the method usually does not know. The area under the curve (AUC) sidesteps the question by aggregating over all thresholds. A perfect classifier fills the unit square and scores 1; the further a curve sags from the corner, the smaller the area. The closer to 1, the better, whatever threshold is eventually chosen.

Choosing that threshold is not a data problem; it comes from the application and the people who own it. A doctor looking at a slightly suspicious exam would rather order more tests. A deforestation alert that dispatches inspectors should not cry wolf. Rio’s SMS rain alerts face both risks: warn too often without rain and people stop believing, fail to warn before a storm and the result can be a catastrophe. The nature of the problem, not the data, should set the threshold.

Conclusion

The lecture’s case against “more of the same” was made in steps. The 0–1 loss has zero gradient almost everywhere, and the squared loss on 0/1 labels produces non-probabilities and punishes confident correct answers. Modelling Y as a Bernoulli replaces both with the logistic loss; making the log-odds linear yields valid probabilities and a convex problem, which Newton’s method solves as a sequence of weighted least-squares fits. Coefficients must be read conditionally, as the student example showed, and the Multinomial extends everything to many classes.

The generative route reached the same destination from the other side: model the features within each class, then apply Bayes’ rule. A shared Gaussian covariance gives LDA and its linear boundaries, separate covariances give QDA, and independent features give naive Bayes. Written as log-odds, all four methods sit in one picture. Logistic regression and LDA are equally flexible, QDA adds interactions, and naive Bayes adds arbitrary per-feature functions, so which one wins depends on whose assumptions match the data. Finally, the 0–1 error rate turned out to be a poor judge of a classifier on its own. The confusion matrix shows which errors are being made, the threshold trades one kind for the other without retraining, and the ROC curve and its area summarise performance across every threshold.

That leaves a practical question the lecture closed on. With LDA, QDA, naive Bayes and logistic regression all on the table, how should we choose among them, and is “use a test set” really the right answer? That is the subject of the next lecture.

Thanks for reading! If you found this useful, consider following for more articles on machine learning and statistical modelling.

Every guarantee of least squares rests on assumptions. This lecture asks what happens when real data breaks them, and what to do about it.

Every image and example was taken from the following sources:

• Lecture 3 — Regressão Linear e Revisão da Semana 1 — Machine Learning 2026 Summer Course at IMPA
• Machine Learning 2026, Prof. Paulo Orenstein

Introduction

These notes cover Lecture 3 of the IMPA Summer Course — Machine Learning 2026 (Paulo Orenstein), one article per lecture.

If you want extra context, the earlier articles in the series follow the course order:

1. Foundations of Statistical Learning: From Rules to Data, and What Makes an Estimate Good

2. Linear Regression: Least Squares, Inference, and How Far a “Simple” Model Stretches

Everything below is self-contained: you do not need to read the other posts to follow this one, as long as you are comfortable with expectations, variances, and enough linear algebra to recognise a matrix inverse and an orthogonal matrix when you see one.

This lecture closes the first week of the course. It opens with a compressed review — the three decisions behind any machine-learning method, empirical risk, the test set, bias and variance, and the whole linear-regression story in one slide — before paying the debt left from the previous class: the five problems that appear when linear regression meets real data (correlated errors, heteroscedasticity, outliers, leverage points, collinearity), each with a diagnostic and a fix. It ends with a question that sounds like an implementation detail and is not: how is β̂ actually computed, if inverting a matrix is a bad idea? The second half of the class was a hands-on Python lab; this article covers the expository part only.

The first week on one slide

Where are we? By the end of Chapter 2 of the course, supervised learning had settled into a small number of moving parts, and it is worth having them in front of us before we start breaking things.

We receive training data — pairs (xᵢ, yᵢ), with features in some space 𝒳 and outputs in 𝒴, usually the real line for regression or a set of categories such as {0, 1} for classification. The goal is a function that predicts well, and “well” is measured by a risk, the expected value of a loss:

In words: among all functions in a class ℱ, pick the one with the smallest average loss. Two things in that line are choices. One is the class ℱ, which we have not committed to yet. The other is the loss L, which says how bad a wrong prediction is. And the data are being treated as random — the object being minimised is an expectation, not a fixed function we could plot.

About the loss there is, as the lecturer puts it, “not much fight.” For regression we use the squared loss, and the moment we do, the target is determined: it is the conditional expectation 𝔼[Y ∣ X = x]. For classification with two classes we use the 0–1 loss, an indicator that is one when we got the class wrong, and the target becomes the most probable class given x. Both are sensible; the first week spent some effort convincing us of that.

Sensible or not, we cannot compute either minimiser, because we have no access to the distribution of Y given X. The way out is to swap the population average for a sample average — the empirical risk — and minimise that instead.

The hat on f̂ is doing honest work: this function was trained on our data and depends on the sample.

That one line is the whole course in miniature, because it exposes the three basic components of any machine-learning method: the class of models ℱ, the loss L, and the optimizer — how we actually find, inside the class, the function that minimises the objective. The strengths and weaknesses of every method we will meet come down to how those three decisions are made, and the third is the subtle one; it will stay subtle for a while. Linear regression picks a very specific class and gets away with a trivial optimizer, because the problem is convex with a closed-form solution. Nearest neighbours picks essentially no restriction on the class and pays for it with an optimizer that is much looser — nothing is being solved exactly.

Judging the choice: test error, and where the error comes from

Once we have a f̂, we need to know whether it is any good. Averaging the squared errors on the training set is a bad idea, and the argument is worth restating: in polynomial regression the training error only goes down as you add higher powers, so a training-error criterion would always pick the highest degree available — taken to its conclusion, it would tell you to include every feature that exists. The fix is a test set, a precious block of data the model never sees, on which the error is computed honestly.

The test error also has a decomposition the lecture calls both universal and intuitive. For a fresh point (x₀, y₀) with y₀ = f(x₀) + ε, add and subtract the average prediction 𝔼[f̂(x₀)] inside the square, expand, and the cross terms vanish:

where the bias is the gap between the average prediction and the truth. The last term is error we commit no matter which method we choose — it is irreducible. The first two are ours. Bias is the error of approximating f with an f̂ that cannot capture its nuances; it comes from a class that is too restricted, linear methods fitting a non-linear truth being the canonical case. Variance is how much the estimate would change if the training data changed; it comes from a method flexible enough to stick to the particular sample it saw. In the polynomial picture: if the truth has degree 10 and we fit a parabola, infinite data will not save us; if we fit degree 50, we interpolate every point and capture noise instead of signal. If we can tell which of the two is happening, we have some hope of improving the model.

That is Chapter 2 in two slides. Chapter 3 applied all of it to one specific class of models, and that is where the lecture turns next.

Linear regression in one slide: β̂ is random, and that is the point

When the class ℱ is the linear functions, choosing f̂ collapses to choosing p + 1 numbers — much easier than choosing an object in an infinite-dimensional space — and the empirical-risk problem has a closed-form answer. Stack the responses into a vector y and the features into a matrix X with n rows and p + 1 columns (the extra column is all ones, for the intercept). Then

There is a lot behind that formula. On the linear-algebra side, it says we are projecting y onto the column space of X. On the statistical side — the part that matters today — betâ is a random object. It depends on y, and y carries a random error, so it makes sense to ask for its mean and variance. If the errors are independent with zero mean and common variance σ²,

This is more useful than it looks. Add the distribution we consider most reasonable for the error — a Gaussian — and we get the distribution of betâ for free:

And once we have the distribution of betâ we can do hypothesis tests and confidence intervals, which the lecturer describes as the most important method behind essentially every scientific discovery of recent years. This is how science decides what is true.

One comma in that story: the tests involve σ², which we do not know. It is estimated by the average squared residual, with a denominator that accounts for the p + 1 parameters already spent:

Besides being reasonable, this is the maximum-likelihood estimator — maximum likelihood being on the course’s prerequisite list, which the lecturer takes as one more reason to revisit that list. With sigmâ² in hand, the test statistics follow: a t test for whether a single coefficient is zero, an F test for whether a whole group of coefficients is zero at once.

The philosophy is always the same, and it will come up repeatedly in the computational part of the course. Assume the null hypothesis. Use a statistic whose distribution under that hypothesis is known. If the value you observe is too extreme for that distribution, then, rather than concluding that unlikely things happen, you conclude that you probably had the wrong view of the world, and you reject. The p-value is the probability, under the null, of seeing something as extreme as what you saw or more; a very small p-value is what triggers the rejection.

A hypothesis test does not prove the alternative. It says: if the null were true, this data would be a surprise — and surprises are evidence that the null is wrong.

That is the machinery. The previous lecture then asked what you actually do with it.

Four practical questions and three extensions

Four questions fall out of this way of seeing the world. The review is brief, because the point of today’s class is what comes after.

Which subset of predictors should I use? You cannot test every subset. The answer was computational: forward selection (start with nothing and add features one by one until the model is good enough), backward selection (start with everything and remove what looks irrelevant), or a mixed strategy. The question will return, because it is really the same question as “how many powers do I include in a polynomial regression?”, and it is a serious problem without a good general answer. In genetics a single test gives you thousands of features per person, and you rarely have thousands of people; linear regression breaks when p > n — intuitively, and mathematically too, since Xᵀ X is then not invertible.

Is any predictor actually significant? Not by checking whether betâⱼ differs from zero — it almost always does. Regress “will it rain tomorrow” on “how many grains of rice I ate today” and a specific sample will give you something like 0.04, because betâ is random and every sample wobbles. The t and F tests are the technology that says whether 0.04 is an effect or just the wobble you would expect around zero if the true coefficient were zero.

Is the fit good enough? The test error is the excellent answer. If you only have training data, training error is not a good metric, but you can weight it against the number of features — the residual standard error does that with its denominator:

The RSE is a mathematical compensation for adding features; the test set is a computational one.

How do I predict, and can I put an interval on it? The prediction is the new feature vector times betâ, and there are formulas for prediction intervals around it. This is the interval people actually want: a doctor does not want an interval on betâ, which they cannot interpret, but “this person has the disease with probability 43%, and with 95% confidence that number is between 40 and 50.”

The extensions were simple in form and large in effect. A binary categorical predictor enters through an indicator — a dummy variable, artificial in the sense that we created it — so the model reads β₀ + β₁ plus noise for observations in class 1 and β₀ plus noise otherwise; with more classes you add one indicator per non-base class. An interaction lets two features talk to each other through their product,

so that the effect of raising X₁ depends on X₂; the hierarchy principle says to keep the main effects when you include the interaction. And non-linearity comes from powers of a feature — polynomial regression — because what is linear in linear regression is the relationship between y and the coefficients, not the features. Given a good basis for your function space, linear regression is exactly the tool for writing y in that basis.

What was left from Chapter 3 is the problems. We made a lot of assumptions along the way and, in the lecturer’s words, the real world is the real world: everything fails. The hope is not that the assumptions are always exactly true but a form of protection: if they fail a little, the methods still basically work; if they fail badly, we have a fix. There are five things to check, the same five as in the course textbook: correlation between the errors, non-constant error variance, outliers, leverage points, and collinearity. We take them in order.

Problem 1: correlated errors

Everything inferential in the last section rested on one line.

Two hypotheses hide in “iid Normal,” and they have very different severities.

Normality is generally not the end of the world. The Gaussian is a very generic distribution in nature: when many essentially random errors accumulate into one phenomenon, the result tends to be Gaussian, which anyone who sat through an engineering statistics course has seen. Still, if it is not precisely true, the inferential procedures — standard errors, confidence intervals, hypothesis tests — become invalid, because each of them used “the distribution is Gaussian, therefore…” somewhere. There are fixes, and the problem sets will get to them.

Independence — the first “i” in iid — may very well be the end of the world. Think about time series in finance, or spatial series in meteorology: data that vary over time, or over space. It would be very strange for the errors there to be independent. If a company is doing very badly today, you do not expect it to recover perfectly tomorrow, nearly go bankrupt the day after, and bounce back again; if things are going well now they will probably keep going well, and if they are going badly they will probably keep going badly. That kills independence between consecutive errors. In space it is the same story: if it is raining at one point, the chance that it is raining just next to it is not low — there is a cloud there — so the error you collect at one location tells you something about the error at the neighbouring one.

In those settings the independence assumption is simply too strong, and the generic remedy is to model the correlation: model enough of it that, once accounted for, it disappears from the errors. Time series has a whole set of methods for this, several of them built on linear regression; spatial statistics has its own; the slides mention Gaussian processes as a specific tool. None of that is covered in this course, though its second half brings methods that handle this kind of structure well — text being the flagship example: if the previous character was “ã”, anyone who has read Portuguese has a good idea of what comes next, and that dependence is exactly what a model like a GPT is built to capture. The point for now is narrower: linear regression can reach these problems, as long as the correlation is modelled rather than ignored.

How do you see correlation in the errors? Plot the residuals against the observation index.

When the residuals go up and down independently, as in the top panel, the iid picture is plausible. When they look like the bottom panel — tendencies that persist, a residual that keeps climbing and then keeps falling — you certainly do not want to apply linear regression naively, at least not for its inferential procedures: t tests, F tests, confidence intervals.

Correlation is about the errors talking to each other. The next problem is about the errors being different sizes.

Problem 2: heteroscedasticity

There was a subtle assumption in the setup that, as the lecturer notes, nobody called him out on. Writing 𝕍[εᵢ] = σ² says the variance of the error is the same for every i. What if it is not? What if, for instance, the errors grow with the features?

His example is a bathroom scale. Weigh people, and the scale has its normal, modest variation. Start weighing much heavier things — a sofa, say — and the variance of the readings may grow with the weight being measured. The variance of the error now depends on the value. That situation has a name that the lecturer clearly enjoys, heteroscedasticity, and it means precisely this: 𝕍[εᵢ] = σᵢ², a variance that changes with i.

Even with independent, Gaussian errors, this can happen, and the diagnostic is a plot of the fitted values ŷᵢ — the prediction for each observation — against the residuals.

Read it as the scale example: when the prediction is small the error is small, and as the predictions get much larger the error grows a lot. Within the instrument’s expected range there is little deviation; at the extremes, the deviation is extreme too. A student asked about the blue curves, and the answer is that they behave like a predictive band on the residuals — capturing, say, the 5th and 95th percentiles, or the 25th and 75th, the lecturer was not sure which — so that even though the mean residual is near zero everywhere, the band makes it visible that the uncertainty is not.

There are two fixes.

The first is to transform the data. Instead of modelling y = Xβ, model the logarithm of y as Xβ. This is still a linear regression — still linear in β — but it changes the nature of what is being predicted. The log barely moves small values and flattens large ones, pushing them down, which stabilises the regression. Here is the same example from the textbook after the transformation.

The second fix is a more parametric modelling of what is going on: weighted linear regression. Ordinary least squares treats every point in the sample as worth the same — the first observation, the second, the tenth. Figure 2 says that is not quite right: the points with high predictions carry a lot of uncertainty and should matter less when we estimate betâ, while the points where we predict 10 and the values really do cluster near 10 are well estimated and deserve priority. Giving each observation a weight formalises that instinct; deriving the weighted estimator is left for the problem sets.

Heteroscedasticity is about the whole cloud of errors changing shape. The next two problems are about individual points misbehaving — first in y, then in x.

Problem 3: outliers

The lecturer’s definition of an outlier is, by his own admission, relatively subjective, like everybody’s: it is a point whose value of y is very improbable. Something with an enormous error, generally for structural reasons. You are trying to predict something that is outside everything you would expect.

The textbook figure shows the same data three ways.

Now the important observation. Points like this do not usually affect the fit. Delete point 20 in your head and ask what happens to the line: it would move from the dashed blue one to the red one, which is basically no move at all. The prediction stays the same. What the outlier does ruin is the evaluation metric. Compute the mean squared error and it explodes, because there is one horrible point in the sum.

So how do you deal with it? This is where the lecture gets, in its own word, polemical.

If you have a deep conviction that the outlier exists because of a structural error — something that certainly should not be there — you remove it. Suppose you are measuring the heights of the students in the room: 1.70 m, 1.80 m, 1.90 m, 1.60 m, and then a 6-metre reading. Surely, at the moment of measuring, someone measured the auditorium instead of a student. That number captures nothing about the height of a person in this course, and it can go.

But notice the word surely. In practice this move is very, very dangerous, because it only makes sense if you can convince everyone else. Imagine publishing your results and writing “this point I did not like, so I threw it away.” It is extremely suspicious: precisely the point where your prediction is terrible is the one you discard, and suddenly your predictions are good. The bar for convincing others that a point can be dropped is very high, and it is not something to do lightly.

If you cannot clear that bar, the outlier is telling you something about the model. Perhaps the model needs to be more complex, or needs other predictors. Or perhaps it is telling you that the error distribution is wrong: under a Gaussian, improbable things do not happen, and you are looking at one. Maybe the errors should be modelled with a heavier-tailed distribution — Laplace, for instance. What happens to the regression when you do that is, again, a problem-set question.

An outlier barely moves the least-squares line, but it wrecks the error metric — and removing it requires convincing the world that it was a mistake, not an inconvenience.

Outliers are unusual in y. Their cousins are unusual in x, and those do move the line.

Problem 4: leverage points

If outliers are unusual values of y, leverage points are unusual values of xᵢ. And they come with a problem the outliers did not have: they definitely affect betâ.

Point 20 does not make much difference to the line. Point 41 does, because the line has to twist to accommodate a point that is far away. Push that point even further out and nudge it slightly, and the entire line would have to move a lot because of one observation. That is, as the lecturer puts it, very uncomfortable.

A second panel makes a subtler version of the same point.

The advantage over outliers — where it was hard to convince anyone a problem even existed — is that leverage has a very reasonable measure. The idea is a derivative: take the prediction ŷᵢ and differentiate it with respect to yᵢ. If I nudge the observed value, does my prediction at that point move a lot? For point 20, no — the prediction there stays essentially put. For point 41 the change is drastic.

To compute the derivative, recall that the fitted values are a linear function of the observed ones, ŷ = Hy, where H is the hat matrix X(XᵀX)⁻¹Xᵀ. Then ŷᵢ is the i-th row of H times y, and differentiating with respect to yᵢ leaves the diagonal entry.

This is the leverage, or influence value, of observation i: how much the fitted value at i responds to its own observed value. It lies between zero and one — proving that is left as an exercise — and it gives you a number to check whether leverage points are showing up in your regression.

Both outliers and leverage points are about single observations. The last problem is about the predictors themselves.

Problem 5: collinearity

Collinearity came up in the previous lecture; here the goal is to add a couple of subtleties. The problem is a linear relationship between predictors.

Take the credit-card example from the course. You want to predict a customer’s end-of-month balance, and you regress it on the credit limit the bank gave them and their credit rating: balance on limit plus rating. Sounds reasonable. What you do not know is that the bank secretly sets the rating as a linear function of the limit — say rating = 80 + 0.06 × limit.

Why is that a problem? Because with those two features the design matrix X no longer has full rank, so Xᵀ X is not invertible — and then, going back to the formula, betâ = (Xᵀ X)⁻¹ Xᵀ y simply does not exist.

The example from last week was the dummy-variable version. You include one indicator for “has a PhD” and another for “does not have a PhD.” The implicit assumption that bites is that X already contains a column of ones, for the intercept. The two indicators sum to exactly that column, so you are back to a linear relationship, and the matrix is again singular. In practice, when you sit down to write code, the software detects this automatically and handles it — in effect, “you are being silly, let me fix that.” It takes some carelessness to land here — an indicator for Brazilian alongside one for non-Brazilian, say.

It is worth being concrete about what actually goes wrong, because it is more than numerical instability in a matrix inversion. Suppose two features are identical, X₁ = X₂, and you ask the model for β₀, β₁, β₂. Predicting with the triple (betâ₀, betâ₁, betâ₂) gives exactly the same predictions as predicting with (betâ₀, betâ₁ + 100, betâ₂ − 100):

and the last term is zero because x₁ = x₂. You are asking the model to do something impossible. You give it two features and ask how much weight to put on each, and it answers that there are infinitely many combinations of weights giving the same prediction, with no way to tell them apart. The parameters are not identifiable. That is the real reason the variance of betâ goes to infinity in this case — and with it, obviously, the hypothesis tests and confidence intervals go out the window.

Perfect collinearity is relatively rare. What actually happens, and is more problematic, is approximate collinearity: two features that are almost linearly related, so that combined linearly they are almost a third one. The model does pick a betâ₁ and a betâ₂, but it starts mixing up the subtle difference between two nearly-identical features with the ordinary statistical error present in any regression. Change the sample a little and the estimates would shift in a different way; attributing effect to one feature or the other becomes very hard. Even in these less extreme cases, the variance of the estimators can get very high.

The genetics example shows why this is not a corner case. With one feature per gene, you would expect several genes to serve the same function — the immune system, say — and to be expressed together. That alone makes the regression suffer from something like collinearity: the estimated variance of betâ comes out enormous, which is the model’s way of saying “you are asking me to distinguish coefficients for features that are all very similar.” Including too many features has this price: beyond being useless for prediction, they can blow up the variance of the estimator. The course will come back to this.

Diagnosing it: from correlation to the VIF

A student asked whether there is a test for which parameters are entangled with which. With just two predictors, plain correlation is enough: in the credit example, every time the rating goes up, the limit goes up too, and the correlation captures exactly that. With p predictors it is more delicate, because pairwise correlations do not reveal multicollinearity — one feature being nearly a linear combination of several others.

The diagnostic is the variance inflation factor. For the j-th coefficient,

where the R² in the denominator is the one you get by regressing Xⱼ on all the other features — using the j-th feature as the response and everything else as predictors. If you are in doubt, you compute this for each feature and look for values that are very large. It is a partial solution, but it does measure the degree to which Xⱼ is explained by the rest.

The formula looks, in the lecturer’s phrase, absolutely crazy, but there is a reason for it. The variance of the j-th coefficient can be decomposed into a factor that depends only on the j-th feature and a factor that depends on all the others.

The first fraction is what the variance of betâⱼ would be on account of Xⱼ alone. The second fraction — the VIF — is how much the other features inflate it. That is what makes the VIF a genuine generalisation of correlation to p variables: it isolates the influence of X₋ⱼ on the variance of betâⱼ. The proof may come later in the course; the slides give the intermediate step, in which the jj entry of the inverse is written in terms of the regression of Xⱼ on the others.

Collinearity is not a numerical nuisance. It means the data cannot tell two explanations apart, and infinite variance is the honest answer to an unanswerable question.

So: linear regression has several problems, and there is a fix or a diagnostic for each one. That is not surprising — Gauss was already looking at this, and there has been time to think of good solutions. One question remains, and it is about the computer rather than the statistics.

Computing β̂ without inverting a matrix

The closed form for the coefficients has an inverse in it, and depending on which path through mathematics you took, that inverse should make you uneasy. Inverting a matrix is a real problem: it tends to be numerically unstable and slow. So the question is whether you have to. Go back to the derivation and the answer is no. There are smarter ways to get betâ than computing an inverse, even on a computer — especially on a computer — and they go through two decompositions from linear algebra, QR and SVD, each useful in different circumstances.

The QR route

Write the design matrix as X = QR, where Q is orthogonal — meaning Qᵀ Q = I — and R is upper triangular. Now go back to the first-order condition, the moment where the inversion happened, and refuse to commit that original sin.

Substitute the decomposition. On the left, Xᵀ X = Rᵀ Qᵀ Q R, and the Qᵀ Q in the middle is the identity and dies, leaving Rᵀ R. On the right, Xᵀ y = Rᵀ Qᵀ y. Both sides start with Rᵀ, which cancels, and the problem collapses to

You could invert R at this point, but you do not even need to, because R is upper triangular, and a triangular system reads off from the bottom up. The last row of R is zero everywhere except its last entry, so multiplying it into betâ kills every term but one, and the last coefficient falls out directly. The row above it has two non-zero entries — one multiplying the coefficient you just found, one multiplying the next unknown — so you isolate that unknown and solve. Keep going up the rows and every coefficient comes out by back substitution, one division at a time. There is no matrix inversion anywhere in the procedure.

The SVD route

Alternatively, write X = U D Vᵀ, with U and V orthogonal and D diagonal. The same trick as before — the inner Uᵀ U is the identity — gives

Now the formula for betâ asks you to invert that, and that turns out not to be hard at all. The inverse of an orthogonal matrix is its transpose, and the inverse of a diagonal matrix is the diagonal of the reciprocals, so the inverse is V D⁻² Vᵀ. Multiply through by Xᵀ y = V D Uᵀ y, and the Vᵀ V that appears in the middle disappears for free.

The only thing that ever gets inverted is D, a diagonal matrix, which is trivial. When you sit down to write code, one of these two routes is what runs underneath any least-squares call — QR or SVD, depending on the sizes involved.

Conclusion

The lecture’s first half was a map of the territory covered so far: every method is a choice of model class, loss and optimizer; empirical risk replaces the expectation we cannot compute; a test set replaces the training error we cannot trust; the test error splits into bias, variance and noise; and for linear regression all of this is explicit, down to the distribution of betâ that makes t tests, F tests and p-values possible.

The second half was the bill for that explicitness. Correlated errors invalidate the inference, and you spot them as persistence in a residual sequence and fix them by modelling the correlation. Heteroscedasticity shows up as a funnel in the residual-versus-fitted plot, and a log transform or a weighted regression closes it. Outliers barely move the line but wreck the error metric, and removing one requires an argument the rest of the world will accept. Leverage points do move the line, and the diagonal of the hat matrix measures exactly how much. Collinearity makes the parameters non-identifiable — infinitely many coefficient vectors, one prediction — and the VIF isolates how much each coefficient’s variance is inflated by the other features. And finally, betâ should never be computed by inverting Xᵀ X: a QR decomposition with back substitution, or an SVD that reduces the inversion to a diagonal matrix, does the job stably.

The recurring lesson is that the assumptions behind a method are not fine print. They are the places where you should expect it to break, and knowing where to look is most of the diagnosis. With the linear model understood — guarantees and failure modes alike — the course is ready for the next family of methods.

Thanks for reading! If you found this useful, consider following for more articles on machine learning and statistical modelling.

Every image and example was taken from the following sources:

• Lecture 2 — Regressão Linear — Machine Learning 2026 Summer Course at IMPA
• Machine Learning 2026, Prof. Paulo Orenstein

Introduction

These notes cover Lecture 2 of the IMPA Summer Course — Machine Learning 2026 (Paulo Orenstein), one article per lecture.

If you want extra context, the earlier article in the series follows the course order:

1. Foundations of Statistical Learning: From Rules to Data, and What Makes an Estimate Good

Everything below is self-contained: you do not need to read the other post to follow this one, as long as you are comfortable with expectations, variances, basic calculus, and enough linear algebra to be unafraid of a matrix inverse.

This lecture instantiates the general supervised-learning setup for the simplest model in the course. We derive the coefficients explicitly, work out their means, variances and distributions, build confidence intervals and hypothesis tests on top of them, rewrite everything in matrix form until a geometric picture appears, prove that this estimator is optimal in a precise sense, and then show that the model bends to accommodate categories, interactions and curves — because “linear” was never a statement about the features.

Where linear regression sits in the learning setup

The general supervised setup has four moving parts: training data, a class of predictive functions, a loss that scores predictions, and an optimizer that picks the best function in the class. Linear regression is what you get when you fill in each part with the simplest reasonable choice.

The training data are pairs, (x₁, y₁) through (xₙ, yₙ), with features in ℝᵖ and responses in ℝ, and we assume n > p throughout — more samples than features. Keep a concrete question in mind while the notation piles up: how many hours of sleep will you get tonight? That is Y; how much magnesium you took, whether you exercised, how old you are — those are the features X.

The predictive class is restricted to functions of the form

In words: a weighted sum of the features plus a constant. That restriction is what gives the method its name.

The loss is the familiar squared error, L(y, f̂(x)) = (y − f̂(x))². If I predict eight hours and you sleep seven, I missed by one, and squaring makes the miss positive either way. The optimizer must then find the coefficients that make the average of that loss as small as possible:

In words: choose the coefficients that minimise the average squared miss on the training sample.

Here is what makes this lecture different from every later one. That optimization problem has a closed-form solution — an explicit formula, no iterative method, no learning rate, no convergence to worry about. Optimization, which will dominate the course once neural networks arrive, is almost invisible here.

That simplicity buys two things. We can understand every mathematical detail of this model, more than we will be able to see once the models get sophisticated enough to hide their internals; and precisely because we understand every piece, the model is easy to extend. Someone who uses this simple method well often goes deeper than someone running a neural network without knowing what is happening inside it.

The plan is four parts: simple regression with one feature, multiple regression with many, extensions of the linear model, and the practical problems that arise when the assumptions we are about to make fail. The last part needs its own lecture — we will assume freely here and pay for it later.

Simple linear regression: finding the line

Start with a single feature. The assumed model is

with β₀ and β₁ unknown numbers we want to estimate and ε the error term. Once we have estimates, prediction for a new point is just ŷ = betâ₀ + betâ₁ x — feed in how much magnesium, read off the hours of sleep.

We choose the coefficients by minimising the residual sum of squares over the training sample:

In words: pick the line whose vertical misses, squared and added up, are as small as possible.

This objective is convex in the two coefficients, so there is no need for anything clever: differentiate, set to zero, solve. Differentiating in the intercept gives a sum that must vanish, and since the intercept does not depend on i, we can isolate it directly:

where the bar denotes the sample average — a piece of notation that will be with us from now on. That expression still contains betâ₁, so it is not yet a solution; we need a second equation. Differentiating in the slope and substituting the intercept we just found gives one, and isolating produces

That is already usable, but it is usually written in a symmetric form, with the deviation from the mean on both sides of the numerator:

The rewrite needs justification. Subtracting ȳ inside the numerator costs nothing because the extra term is exactly zero: ȳ comes out of the sum, and what is left is the sum of xᵢ − x̄, which vanishes by the definition of x̄. The same argument, with the roles of xᵢ and x̄ swapped, kills the corresponding term in the denominator. Small identities like this keep reappearing, so it is worth doing on paper once — and since slides make all of this look like it follows from nothing, this is the one lecture where redoing the algebra yourself is genuinely within reach.

> Both formulas are computable from the training sample alone. Nothing in them refers to a quantity we do not have.

Now put it to work on the advertising data, regressing product sales on the TV budget. The formulas give betâ₀ = 7.03 and betâ₁ = 0.0475, so the fitted model is

The slope has a direct reading: every 100 dollars put into TV advertising is associated with an estimated increase of 100 × 0.0475 = 4.75 in sales.

The fit looks reasonable and the increasing relationship is obvious. But the model pays a price at the extremes, and the variability of the points grows with the budget — hold on to that, it is one of the assumptions we will have to revisit.

What is random here, and what it costs us

Everything so far treated the estimates as numbers. They are not. Both are computed from the training sample, and a different sample would give different values — they are random variables, so what are their mean, variance and distribution?

Fixing the features, the only source of randomness is the error term, with 𝔼[εᵢ] = 0 and 𝕍[εᵢ] = σ². Since the errors average out to zero, 𝔼[yᵢ] = β₀ + β₁ xᵢ, and because expectation is linear and the features are fixed, it passes straight through the estimator’s formula. The β₀ terms cancel in the numerator, a factor of β₁ comes out, and numerator and denominator become the same sum:

In words: on average, the estimators hit the true values. They are unbiased.

Being right on average is not the whole story: an estimator that is right on average but swings wildly can be far worse than one that barely moves. So we need the variance too, and variance being a quadratic operator rather than a linear one makes the computation heavier. Write the slope estimator as a linear combination of the responses, with weights

and use that each response has variance σ²:

In words: the more spread out your feature values, the more precisely you can pin down the slope.

The intercept takes one more step, since it mixes two random quantities. Expanding the variance of a difference gives a term for ȳ, a term for the slope, and a covariance between them. The first is σ² / n; the third is zero, again because the deviations from the mean sum to zero. What is left is

These formulas are exact. They are also useless as they stand, because both contain σ² — the variance of the error — and we do not know it. It has been a symbol this whole time. Is it 10? Is it 0.4? Nobody knows, which means nobody can actually compute the standard error of a coefficient from a sample.

The way out is the obvious one, and it is where the argument starts folding back on itself: estimate it.

The intuition: if the estimates are close to the true coefficients, the residual yᵢ − ŷᵢ is close to the error εᵢ itself. And since the errors have mean zero, their variance is the expectation of their square — approximated here by the sample average of the squared residuals. Why the divisor is n − 2 rather than n is genuinely not obvious, and the lecture leaves it as a debt to be settled once we have the matrix machinery.

The distribution, and why we need one more assumption

Mean and variance are not enough to build a confidence interval; for that we need the distribution. And that requires an assumption we have not made yet: that the errors are independent and normally distributed,

This is a genuine addition: before, we only assumed the errors had mean zero and variance σ²; now we are pinning down their shape.

The payoff is immediate. A linear transformation of a normal variable is still normal — the mean and variance may change, the shape does not — and both estimators are linear transformations of the responses, hence of the errors. So both are normal, and a normal is fully determined by its mean and variance, which we already have. Standardising the slope,

and since a standard normal satisfies ℙ[|Z| ≤ 2] ≈ 0.95, isolating β₁ inside that inequality gives a 95% interval.

Except the interval contains the standard error, which contains σ², which we do not know. We have to substitute the estimate — and substituting changes the distribution. To describe how, two definitions.

If Z₁, …, Zₙ are independent standard normals, the sum of their squares follows a chi-squared distribution with n degrees of freedom:

And if a standard normal Z and a chi-squared K with n degrees of freedom are independent, the ratio below follows a t distribution with n degrees of freedom:

Neither is mysterious in practice: both are distributions a computer can plot, and any probability you want from them is a lookup.

The fact that ties this together — stated here, proved later in full generality — is

Now divide the standardised slope by that quantity, appropriately normalised. The numerator is a standard normal, the denominator is the square root of a chi-squared over its degrees of freedom, and the two are independent (also proved later). By definition, the ratio is a t:

Notice the two hats stacked in the denominator. The first says betâ₁ is an estimator, so it has a standard error; the second says that standard error itself had to be estimated.

For large n the t distribution is very close to the normal, and ℙ[|t| ≤ 2] ≈ 0.95 still holds. So the interval survives the substitution:

and every quantity in it is computable from the sample. For sales on TV, where betâ₁ = 0.0475, the 95% confidence interval is [0.042, 0.053]. Instead of “your friend will sleep 7:15”, we can say “with 95% confidence, between these two numbers” — a considerably more honest statement.

One caveat, the one everyone gets wrong the first time. It is tempting to read that interval as “β₁ has a 95% chance of being in here”. It does not. β₁ is not random — it is a fixed number we happen not to know. What is random is the interval: draw many samples, compute the interval each time, and about 95% of them will contain the true β₁. The coefficient stays put; the interval moves around it.

Hypothesis tests and the p-value

There is a useful observation buried in that interval. If it does not contain zero, the data are telling us something: β₁ is very likely inside, zero is outside, so β₁ is probably not zero. And β₁ = 0 is precisely the statement that X has no effect on Y.

That reframing is how science usually asks its questions. Formally,

Does the vaccine work? Does the Higgs boson exist? Does magnesium affect sleep? Every one of these, once you have a model, becomes a question about whether a coefficient is zero.

Under the null hypothesis, β₁ = 0, so the statistic

should follow tₙ₋₂ — it is the same statistic as before with the hypothesised value plugged in. Compute it on the observed sample, call the result the observed t, and ask: if we really lived in the world where H₀ holds, how probable would it be to see something at least this extreme? That probability is the p-value:

The logic is a proof by discomfort. Assume we live in the null world; compute how surprising our data would be there. If the answer is “extremely surprising” — say, less than a 5% chance — and yet here the data are, the sensible conclusion is that we probably do not live in that world. So we reject it.

The p-value is a controversial object these days, and deservedly scrutinised, but it remains useful. Applied to the advertising data:

Both statistics are large and both p-values are tiny, so it is implausible that either coefficient is zero: there is a real relationship between TV spend and sales. Which the scatterplot already suggested — and that is exactly the point of the machinery. Your eyes work beautifully in two dimensions and not at all in five, let alone a billion; these formulas are what your eyes are doing, written down.

One honest limitation, raised in class: rejecting H₀ assumes the functional form is reasonable. The whole test is conditional on the relationship being linear, and if it is not, the test is answering a question about a model that does not describe the data. We will get to that.

Measuring the fit: RSE and R²

Before adding more features, we need a way to say how good a fit is. Two measures show up.

The first is the residual standard error, essentially a standardisation of the mean squared error:

It is not used heavily in practice, but it is standard vocabulary and it will come back shortly with a real role to play.

The second is dimensionless and much more common — the coefficient of determination, R², the fraction of the variance of Y that the model explains:

The comparison is against a deliberately silly model: one that ignores X entirely and always predicts the sample mean. Asked how long someone will sleep, this model shrugs and says “in my sample, humans sleep about eight hours”. R² asks how much better than that you did. Predict every point perfectly and R² = 1; do exactly as well as the silly model and R² = 0. So 0 ≤ R² ≤ 1, larger being better.

Multiple regression, matrices, and a picture

One feature is limiting. In the advertising example it makes far more sense to use TV, radio and newspaper together than to pick one. The model becomes

Where the first half of this lecture leaned on probability, this half leans on linear algebra — writing indices by hand for an arbitrary number of features is not a life anyone wants.

Stack the responses into a vector y of length n and the features into a matrix X with n rows and p+1 columns, whose first column is all ones. Rows are observations — your friends — and columns are features: milligrams of magnesium, hours of exercise, age. The leading column of ones absorbs the intercept, and the n separate equations collapse into one:

The objective collapses too — squaring each residual and summing is exactly multiplying the residual vector by its own transpose:

Differentiating with respect to a vector makes people nervous the first time, but two rules cover it: the derivative of aᵀβ is a, and that of βᵀAβ is 2Aβ for symmetric A. Applying them,

Assume X has full rank. Then XᵀX is positive definite, the problem is convex, and setting the gradient to zero gives the first-order condition Xᵀ(y − Xbetâ) = 0, which we solve to get

Strange-looking the first time, and worth getting used to — it is the formula the rest of the summer is built on, and everything we did by hand in the one-feature case is a special case of it.

Predictions follow immediately. Substituting,

H is the hat matrix, for the obvious reason: hand it the observed responses and it hands back the fitted ones. It puts the hat on.

That is a lot of structure for one matrix, and structure is never an accident. H is symmetric, straight from its form, and idempotent — applying it twice does nothing new, since the inner (XᵀX)⁻¹XᵀX collapses to the identity:

Symmetric and idempotent together are the defining properties of an orthogonal projection — here, onto the column space of X.

The geometry makes this concrete. Rewrite the first-order condition using Xbetâ = Hy: it says Xᵀvarepsilon̂ = 0, where varepsilon̂ = y − ŷ is the residual vector. In other words the residuals are orthogonal to every column of X, hence to the whole subspace they span.

So here is what fitting a linear model actually is. The column space of X holds every prediction the model is capable of producing. The response vector lives in ℝⁿ and, since that subspace has dimension p+1 with n > p+1, y generally sits outside it. Linear regression finds the closest point inside: the orthogonal projection. The residual, being the leftover, is perpendicular to everything the model could have predicted.

Least squares is not an optimization problem in disguise. It is a projection, and the normal equations just say that the leftover points straight out of the model’s reach.

The lecture marks this geometry as optional material — not exam content, but the thing that makes the formula stop looking arbitrary. Two assumptions are quietly doing work: full rank of X, and the convexity that follows from it. Convexity, note, is in the coefficients, not in the features. Both come back later.

Gauss–Markov: why this estimator and not another

Nothing so far explains why we should prefer this estimator over the many others we could write down. The Gauss–Markov theorem does.

Among all linear unbiased estimators of β, the least-squares estimator has the smallest variance. Precisely: for any vector c and any estimator of the form betã = Ay with 𝔼[betã] = β,

In words: no linear unbiased estimator beats least squares on variance, in any direction you care to look. This is the property summarised by calling betâ the best linear unbiased estimator, or BLUE.

The proof is short enough to sketch. Write any competing linear estimator as betã = ((XᵀX)⁻¹Xᵀ + D)y for some matrix D — that covers everything linear. Demanding unbiasedness for every possible β forces DX = 0. Expanding the variance with 𝕍[y] = σ² I, the cross term carries a factor of XᵀDᵀ = (DX)ᵀ = 0 and vanishes, leaving the competitor’s variance equal to ours plus σ² |Dᵀc|² — a squared norm, which cannot be negative.

It pays off when connected to the decomposition of prediction error into squared bias plus variance. Our estimator has zero bias and, among unbiased linear estimators, the smallest variance — so it also has the smallest mean squared error in that class. Within the linear unbiased world, this is the right answer, not merely a convention.

Inference when there are many features

With the matrix formulation in hand, every inference result from the one-feature case generalises, and the debts we left behind get paid.

Unbiasedness first. Substituting the model into the estimator,

and the second term has expectation zero, so 𝔼[betâ] = β. The variance follows from 𝕍[ε] = σ² I:

This mirrors the simple case — σ² over a sum of squared deviations becomes σ² times an inverse matrix, with x̄ absorbed into the column of ones. Assuming normal errors, betâ is a linear transformation of a normal vector and is therefore jointly normal:

The error variance is estimated the same way as before, with the degrees of freedom generalised:

Setting p = 1 recovers the n − 2 from the simple case.

For a single coefficient, write betâⱼ = eⱼᵀbetâ, where eⱼ has a one in position j and zeros elsewhere. Sandwiching the covariance matrix between two copies of eⱼ just picks out a diagonal entry:

Paying the chi-squared debt

Now the result we assumed on faith earlier. Start from the residual vector and use HX = X, which kills the systematic part:

Since I − H is symmetric and idempotent for the same reasons H is, the sum of squared residuals is a quadratic form in the errors alone: (n − (p+1))sigmâ² = εᵀ(I − H)ε.

Take a spectral decomposition I − H = QΛ Qᵀ with Q orthogonal. Idempotence forces every eigenvalue to satisfy λ² = λ, so each is zero or one. Rotating the errors by an orthogonal matrix leaves their distribution alone — the spherical normal is rotation-invariant — so the rotated errors are still independent normals with variance σ². The quadratic form is therefore a sum of squared standard normals, with as many terms as there are eigenvalues equal to one:

And the trace is easy. The trace of a projection matrix is the dimension of the subspace it projects onto, so tr(H) = p+1, giving tr(I − H) = n − (p+1). Hence

The n-2 version we used earlier is the case p = 1, and the mysterious divisor was never mysterious: it counts the dimensions the projection leaves over. One loose end — the t distribution requires numerator and denominator to be independent, and betâ and varepsilon̂ are, because the covariance between them vanishes. (Left as homework in class; the computation is two lines.)

Testing one coefficient, and testing several

The t test now carries over unchanged apart from the degrees of freedom:

But having many features raises a question that could not be asked before. With one feature, the only thing to test was whether that coefficient was zero. Now we might want to know whether several are jointly zero — whether anything at all among magnesium, exercise and diet affects your sleep. No single t test answers that, and running one per coefficient is not a substitute.

The tool is the F test. To test H₀: βₖ₊₁ = ⋯ = βₚ = 0, compare the residual sum of squares of the restricted model — the one fitted with those coefficients forced to zero, giving predictions ŷᵢ⁽ᵏ⁾ — against that of the full model:

Two special cases matter. Taking k = p-1 tests a single coefficient, and the F test then turns out to be equivalent to the t test — a non-trivial but very doable exercise. Taking k = 0 tests whether all of the coefficients are zero, which is the model’s global significance test.

The proof reuses everything we just built. The F distribution is a ratio of independent chi-squareds, each divided by its degrees of freedom. Partition the design matrix into the columns being tested and the rest, and let H₋ₖ be the projection onto the restricted model’s column space. The full model’s residual sum is

and an analogous computation shows that the difference of residual sums is

which is, geometrically, the squared length of the extra projection gained by adding the extra columns.

Independence of the two makes the ratio an F, and it comes from (I-H)(H − H₋ₖ) = 0: projecting onto the big subspace and then onto a smaller one inside it is the same as projecting straight onto the smaller one. Zero cross-covariance plus joint normality gives independence, and the ratio is F by definition.

Reading a regression table without fooling yourself

Time to actually fit the advertising data with all three channels.

Newspaper is not significant, which may be surprising. Here is a clue: the correlation between newspaper and radio spend is 0.35. The two plausibly reach a similar audience, so once money is already flowing into radio, adding newspaper does not buy much — and radio, being strongly significant, absorbs the shared effect.

This is where the interpretation of a multiple-regression coefficient matters. βⱼ measures the average effect of Xⱼ on Y holding all the other predictors fixed — not the number you would get by regressing Y on Xⱼ alone. “Given that I am already spending heavily on radio, does a bit more newspaper help?” is the question the table answers, and the answer is: not much. That is a different statement from “newspaper advertising is irrelevant”.

And the largest trap is one significance cannot protect you from. Take a beach in Australia, regress shark attacks on ice-cream sales, and the coefficient comes out beautifully significant. Ice cream does not cause shark attacks: on hot days more people buy ice cream and more people go into the water, which raises the chance of an encounter. Sell no ice cream and the beach is empty.

A statistically significant coefficient is evidence of association. It is never, by itself, evidence of causation.

The same caution applies to the running example: even a highly significant magnesium coefficient would not show that magnesium causes better sleep, only that the two move together in the data.

Four questions you actually ask in practice

Four questions organise the practical use of multiple regression: is at least one predictor significant? Is there a subset worth keeping? How good is the fit? And how accurate is a prediction for a new point?

Significance. Use the F test with k = 0 to test whether all coefficients are zero at once — under the normality assumption on the errors. To test a single coefficient, put that variable last and use k = p-1, which is the t test again.

Variable selection. With p predictors there are 2ᵖ subsets, so exhaustive search dies quickly. Three practical strategies:

- Forward selection — start from the null model and add, at each step, the variable that most reduces the residual sum of squares: fit every one-feature model (magnesium alone, exercise alone, age alone), keep the best, then try adding each remaining feature to it.
- Backward selection — start from the full model and drop, at each step, the variable with the largest p-value.
- Mixed selection — like forward selection, but a variable already included can be dropped if its p-value crosses a threshold later on.

All three are greedy, and none is guaranteed to find the best subset — a reasonable path, not an optimal one. Better strategies come later in the course.

Fit. There is a trap here that echoes the training-error story: the residual sum of squares always (weakly) decreases when you add a predictor, even a useless one, so judging models by it always picks the fullest one. The residual standard error corrects for this by dividing by the remaining degrees of freedom,

so adding a barely useful variable lowers the numerator but also lowers n − (p+1), and the RSE can increase. Complexity gets penalised rather than rewarded.

Beyond any single number, plot the residuals — against the fitted values, or against each predictor. In the advertising fit they show a distinct belly: over-predicting in the middle of the range, under-predicting at both ends. That shape is the model telling you it is missing a curvature.

Prediction intervals. For a new point, the prediction is ŷₙ₊₁ = xₙ₊₁ᵀbetâ. The prediction error has mean zero and variance

where the first term is the new observation’s own noise and the second is our uncertainty about the coefficients. Assuming normality, the standardised prediction error follows tₙ₋₍ₚ₊₁₎, so with roughly 95% probability

This is a different animal from a confidence interval. There, the interval covered a fixed unknown parameter; here it covers a future observation, which stays random even if you knew the coefficients exactly, because it carries its own error term. “This person will sleep 9 hours, and with 95% probability between 8.5 and 9.5” — that is a predictive statement, and linear regression hands it to you almost for free. Not every method does, and if you are trying to learn something about the world, knowing how much to trust a prediction matters more than the prediction itself. (A note on the word, since it collides with other usage: throughout this course, inference means statistical inference — not what the term means when people talk about running a trained model.)

Extension one: predictors that are not numbers

Everything so far assumed the features are real numbers. Plenty are not. What sport do you play — tennis, judo, swimming, none? That is not a number and cannot go into X as it stands.

Take a concrete setup: predicting credit-card balance from six quantitative variables and four qualitative ones — own (do you own your home?), student, status (marital status), and region (north, west or south zone).

For a two-level variable, define a dummy: xᵢ = 1 if person i owns a home and 0 otherwise. A property of a person has become a number that fits in the matrix, and the interpretation is neat. The model yᵢ = β₀ + β₁ xᵢ + εᵢ now reads as two cases: β₀ + β₁ plus noise for owners, β₀ plus noise for everyone else. So β₁ is exactly the average impact of home ownership on the balance.

On the book’s data, betâ₀ = 509.80 is the average balance for non-owners, and owners come in at 509.80 + 19.73 = 529.53. Tidy — but the p-value on that difference is not small, so there is no evidence the gap is real.

More than two categories needs more than one dummy: one fewer than the number of categories. For region, define xᵢ₁ = 1 for the south zone and xᵢ₂ = 1 for the west zone, and the model reads as three cases — β₀ + β₁ for south, β₀ + β₂ for west, plain β₀ for north.

The category left without a dummy is the baseline. Here that is the north zone: β₀ is its reference level, and the other coefficients measure each category’s differential relative to it. On the book’s data, betâ₀ = 531.00, with -18.69 for south and -12.50 for west, both relative to north, and neither individually significant.

But asking whether region matters at all is not two t tests. It is the joint hypothesis that both coefficients are zero — an F test — and here it returns a p-value of 0.96. Nothing to reject; no evidence that region affects the balance.

One question remains: why leave a category out? Adding features unnecessarily is a fair objection, but the real problem is worse. If everyone belongs to exactly one category, the dummies for all categories always sum to one — which is exactly the intercept column. The matrix loses full rank, XᵀX stops being invertible, and the formula we built everything on has nothing to compute. This is collinearity, and the exact version here is the easy case to avoid. Near-collinearity — features almost, but not quite, linear combinations of others — causes trouble too, and waits in the next lecture.

Extension two: interactions and curves

The additive model says the effect of each feature is the same regardless of the others. In the advertising case, sales = β₀ + β₁ × TV + β₂ × radio + ε implies that spending one more unit on TV always changes sales by β₁, no matter what radio is doing.

The data disagree. When the budget is split between TV and radio, the additive model under-predicts sales; when it is concentrated in one channel, it over-predicts. That is the signature of an interaction: the return on each channel depends on the other — plausible enough, since channels reaching overlapping audiences should have diminishing joint returns, while one reaching a fresh audience is worth more.

The fix is to add a multiplicative term as a new feature:

And the improvement is not merely statistical: R² climbs from 89.7% to 96.8%. The interaction buys real predictive power.

Interpreting it takes a small derivative. Since the partial derivative of sales with respect to TV is now β₁ + β₃ × radio, an extra 1000 dollars in TV is worth 19 + 1.1 × radio in sales, and symmetrically 1000 dollars more in radio is worth 29 + 1.1 × TV. The return on each channel grows with what you already spend on the other.

A rule comes with interactions, the hierarchy principle: whenever you include a product of two predictors, include the individual predictors too, even if their own coefficients are not significant. Without them the interaction coefficient means something else — a model with only student × own and no separate terms would have that coefficient comparing student homeowners against everyone else combined, rather than isolating the interaction.

Interactions also work between a quantitative and a qualitative predictor. Predict balance from income and the student dummy, first without the product term student × income and then with it.

Algebraically, the non-student line is betâ₀ + betâ₁ × income, and without the interaction the student line is the same slope shifted by betâ₂. With it, the student slope becomes betâ₁ + betâ₃ — the effect of income can now be weaker or stronger depending on the group. A far more flexible model for the price of one extra column.

Linear in the parameters, not in the features

The last extension goes after the assumption that has been nagging since the start: why should anything be linear? Nothing forces us to stay linear in the features. We can model

by treating powers of X as additional features. Same machinery, bigger matrix.

The Auto data make the case. Plot fuel economy against engine horsepower and the linear fit is visibly wrong: below the data at low horsepower, above it in the middle, below again at the high end. Adding a quadratic term,

bends the fit to follow the shape of the cloud.

That degree-5 curve is a warning, though: it fits the training data more closely and starts to wiggle, the visual signature of overfitting, which takes us straight back to the trade-off between bias and variance.

The general point is bigger than polynomials, and it is the one worth carrying out of this lecture:

“Linear” in linear regression means linear in the parameters, not in the features. Any transformation of the features — powers, logarithms, products — keeps the problem a linear regression with a larger design matrix.

Which has a striking consequence: expand a feature in a rich enough basis of powers and you can approximate any reasonable function arbitrarily well. The model is far more expressive than its name lets on.

Conclusion

The through-line of this lecture is that a model simple enough to understand completely is a model you can push. We derived the coefficients in closed form, found their means and variances, added a normality assumption and got distributions, confidence intervals and hypothesis tests out of it. We rewrote everything in matrices, watched a projection appear, and proved via Gauss–Markov that this estimator is the best among linear unbiased ones. Then we paid our chi-squared debt using nothing but the geometry of that projection, built the F test on top of it, and showed that categories, interactions and curvature all fit inside the same framework.

The ideas that survive are worth naming: least squares is a projection; degrees of freedom count the dimensions it leaves over; a t test asks about one coefficient, an F test about several; significance is association, never causation; and linearity is a claim about parameters, not features.

What is left over is the bill. We assumed the design matrix has full rank and that the errors are independent, normal, and share a common variance — and we saw hints, in the widening spread of the first scatterplot and the belly in the residual plot, that reality does not always cooperate. Those assumptions and what to do when they fail are the subject of the next lecture, along with the problems that come with them: correlated errors, heteroscedasticity, leverage points, and collinearity.

Thanks for reading! If you found this useful, consider following for more articles on machine learning and statistical modelling.

Every image and example was taken from the following sources:

• Lecture 1 — Introdução ao Curso / Noções gerais de aprendizado — IMPA Summer Course
• Machine Learning 2026, Prof. Paulo Orenstein

Introduction

These notes cover Lecture 1 of the IMPA (Institute for Pure and Applied Mathematics) Summer Course — Machine Learning 2026 (Paulo Orenstein), one article per lecture.

This is the first article in the series, so there is nothing to read before it. Everything below is self-contained: the only background assumed is comfort with expectations, variance, and basic calculus.

The lecture runs one long argument from beginning to end. It starts with a paradigm shift — the move from writing rules to learning them — then makes that goal mathematically precise as estimating a target function *f*. From there it asks what we do with an estimate, how the methods that produce one split into two families, what we give up when we make them more flexible, how we tell a good estimate from a bad one, and what exactly the error we measure is made of. It closes by showing that the whole structure survives when the thing we want to predict is a category instead of a number.

From rules to data: the paradigm shift

For most of computing history, the arrangement was simple. A programmer writes rules. The program receives data. Applying the rules to the data produces answers.

You are relying on this right now. When you click an icon on your desktop and the program opens, it opens because someone wrote the rule connecting that click — the data — to that action — the answer. The same is true of a Breakout-style game, where every behaviour is spelled out by hand: if the ball collides with a brick, the brick disappears and the ball reverses direction.

This works beautifully until you hit a task where you know how to do something but cannot say how you do it. Driving a car, for instance. When do you turn? How do you react to a pigeon in the road, a pedestrian stepping off the curb, rain on the windshield? You do all of this competently and you cannot write it down as a finite list of conditions. Or take an easier-sounding task: deciding whether a photo contains a cat. Two ears? Plenty of things have two ears, and some cats in some photos show one. Whiskers? Sometimes invisible. Every feature you propose comes with exceptions, and the exceptions come with exceptions.

Machine learning inverts the arrangement. Instead of supplying rules and data to get answers, we supply data and answers, and the program’s job is to produce the rules.

Stated formally, and this is the sentence the entire course unpacks: we have data X, we have answers Y, and we want to learn a function f such that

In words: find the rule that turns our inputs into our answers.

That single reformulation is what put machine learning everywhere. Consider satellite imagery, where each pixel carries four spectral bands. Those bands are X; the predicted land use for that pixel is Y. Nobody hand-codes the spectral signature of farmland versus water — the labels teach it.

Recommendation works the same way. Deciding what film to suggest next used to mean asking someone with taste; now a learning algorithm reads your history and decides.

The same framing also covers questions that are not really about prediction at all — relating demographic factors like age, calendar year, and education level to salary, for instance, where what we want is the shape of the relationship rather than a forecast.

And it stretches all the way to generation. Text-to-image systems take a description and return a picture that matches it.

Chatbots go further still, explaining why a meme about GPUs and VRAM is funny, even though what the algorithm actually receives is a grid of zeros and ones with no visual understanding attached to them.

Prediction, recommendation, inference, generation, conversation — the applications look wildly different, and every one of them is a version of learning f from X and Y.

What is f? The target as a conditional expectation

Saying “learn f such that f(X) ≈ Y” is a goal, not a definition. If we want to do mathematics with it, we need to know precisely what f is.

Take a concrete problem the lecture returns to repeatedly: predicting product sales from advertising spend across three channels. Here Y is sales — the response, or target. The predictors are X₁ for TV budget, X₂ for radio, X₃ for newspaper; together they form the feature vector X.

The working hypothesis about how the data are generated is

where f is an unknown function and ε is independent noise with mean zero. Sales are driven by advertising through some fixed but unknown mechanism, plus randomness we cannot account for. The goal is to discover f.

A picture of exactly that decomposition, in one dimension, makes the pieces concrete.

Now take conditional expectations on both sides of the model, given X:

Two facts do all the work. Once we condition on X, the quantity f(X) is no longer random — it is a fixed number — so its conditional expectation is itself. And the noise has mean zero, so its conditional expectation vanishes. What remains is a clean characterisation:

In words: the function we are chasing is the average value of the response among all cases sharing that particular input. With three predictors it reads f(x₁, x₂, x₃) = 𝔼[Y | X₁ = x₁, X₂ = x₂, X₃ = x₃] — the mean sales among all campaigns with exactly that budget split.

The target function is not an arbitrary curve through the data. It is the conditional mean of the response given the inputs.

There is a second route to the same destination, and it is worth taking because it does not require the additive-noise hypothesis at all. Suppose we make no assumption about how Y was generated, and simply decide that a good predictor is one with small squared error. Then among all possible functions g, the one minimising expected squared error is

The argument is a short calculation. Expanding the square gives

which is a quadratic in the number g(x). Differentiate with respect to g(x), set the derivative to zero, and out comes g*(x) = 𝔼[Y | X = x].

So the conditional mean is the best predictor under squared error whether you arrive at it by assuming Y = f(X) + ε or by simply choosing squared error as your criterion of success. Two different starting points, one answer.

In one dimension your eye can trace that curve straight off the scatter plot. The formalism earns its keep the moment there are more variables than a plot can hold — and that is the ordinary case, not the exception.

And now the catch that shapes everything after it. We do not know the distribution of Y given X. If we did, we would compute 𝔼[Y | X = x] and be finished. The target function is perfectly well defined and completely inaccessible; every method in this course is a strategy for estimating it from a finite sample.

Prediction and inference: two reasons to want f

Before building estimators, it is worth asking what we intend to do with one. Write f̂ for our estimate of f. There are two distinct motivations, and they pull in different directions.

The first is prediction. Given a new input X, we compute ŷ = f̂(X) and use it. Spend 100 on TV, 20 on radio, nothing on newspaper — how much will we sell? Here f̂ is a machine for producing numbers, and nobody needs to look inside it.

The second is inference. Here f̂ is not a machine but an object of study, and the questions are about its structure. Which features actually matter? Is there evidence that f is linear? Is money spent on TV more profitable than money spent on radio? Does newspaper advertising do anything at all? These are second-order questions relative to prediction — they compare and explain effects instead of producing a forecast.

One clarification, since the word is overloaded: “inference” here means asking second-order questions about a fitted f̂. It is not the LLM sense of the word — running a trained model forward to produce output — which the course discusses separately later.

The distinction matters because methods are not neutral between these goals. A method that predicts superbly may tell you nothing about which features drove the prediction, and that trade-off has a section of its own below.

Two families of methods: parametric and non-parametric

Whatever the motivation, f̂ has to be learned from a finite sample: training data {(x₁, y₁), …, (xₙ, yₙ)}, where each xᵢ = (xᵢ₁, …, xᵢₚ)ᵀ collects p features for observation i.

The lecture offers a useful analogy for what this learning looks like. Nobody teaches a child what a tree is by supplying a definition — no rule set about bark, branching structure, and minimum height. You point at things and say “tree”, “not a tree”, over and over, and eventually the child gets it. Learning algorithms are in the same position: they see labelled examples and produce a rule.

There are two broad strategies for doing so. Parametric methods make a strong hypothesis about the form f takes. Non-parametric methods make no explicit hypothesis about its form at all. If the parametric hypothesis is approximately valid, parametric methods tend to win. With a lot of data, non-parametric methods can win instead, since no structure has to be imposed and the data alone decide the shape of the function.

The parametric route: linear models

The canonical parametric assumption is linearity:

with unknown coefficients β₀, …, βₚ ∈ ℝ. Look at what this buys us. Estimating an arbitrary function means searching an infinite-dimensional space; estimating this one means finding p+1 numbers. Prediction is then

Linear models are easy to interpret and easy to extend. A straight-line fit f̂(X) = betâ₀ + betâ₁ X may already be reasonable; adding a betâ₂ X² term can fit better still.

But there is an obvious temptation here, and it is a trap. If X² helps, why not X³, X⁴, X⁵? Keep going and the curve will eventually pass exactly through every training point — and that, as we will see shortly, is precisely what we do not want.

The limits of the linear hypothesis show up clearly in a simulated example where the truth is known. Take a surface generating income from two variables, `income = f(education, seniority) + ε`.

Fit a linear regression plane to data from this surface, and the fit gets the broad story right: more education raises income, and seniority matters less. What it cannot do is bend.

The non-parametric route: local averaging

The other strategy starts from the characterisation we derived earlier. If f(x) = 𝔼[Y | X = x] is a population average, then approximate it with a sample average — over the observations near x:

where N(x) is a neighbourhood of x. To predict at a point, look at the training points nearby and average their responses. No functional form is assumed anywhere.

This works well when the number of features p is small and the number of observations n is large — say p ≤ 4 and n ≥ 1000. When p is large it performs poorly, a failure we will return to in later lectures.

More sophisticated non-parametric methods follow the same philosophy. A thin-plate spline imposes no pre-specified form for f either, and because the data decide the shape, it can capture what the plane could not.

Freedom has a price, though, and here it is the neighbourhood size. Make the neighbourhood too wide and every prediction approaches the global average: the fitted curve flattens out and stops responding to x at all. Make it too narrow and the curve clings to individual data points, wobbling up and down, until in the extreme it interpolates every point perfectly.

That last sentence should sound familiar. It is the same failure as the runaway polynomial, reached from the opposite direction. Both families can destroy themselves with excess flexibility, which is exactly the trade-off the lecture turns to next.

Flexibility, interpretability, and what we give up

The lecture names three trade-offs explicitly, and they are related but not identical.

The first is fair fit versus under-fit or over-fit: how do we know when a fit is good? The second is prediction accuracy versus interpretability: linear models are simple to interpret, thin-plate splines are not. The third is parsimony versus black box: a preference for simple models with few variables over complex models with many.

The second of these is not a minor technical footnote. Modern large language models are extraordinarily capable and correspondingly hard to interpret or control precisely — even the people building such systems may be unable to guarantee that a specific behaviour will change. “I can’t”, as the honest answer goes.

Compare that to linear regression, where interpretation is nearly mechanical. Want to know what happens to predicted salary if years of education increase slightly? The coefficient β₁ is the answer, directly and quantitatively. Ask the same question of a spline or a nearest-neighbours estimate and there is no coefficient to read — you have to probe the fitted function and reason about what it does locally.

Parsimony has a similar flavour. In a genetic panel used to predict cancer risk, thousands of variables are measured and only some are believed to be genuinely relevant. Finding which ones matter is already hard; it gets harder as the model gets more complex. A model with few variables is not merely more elegant — it is more answerable.

Putting methods on a plane with flexibility running left to right and interpretability running bottom to top produces a fairly consistent ordering.

The relationship is roughly inverse, and it explains the shape of the course itself. The first half stays in the linear-model world — with extensions like splines and generalized linear models — precisely because that world remains interpretable. The second half moves into the high-flexibility, state-of-the-art methods that are much harder to understand from the inside.

Which leaves the first trade-off still open: how do we tell a good fit from a bad one?

Measuring accuracy: training MSE, test MSE, and overfitting

For a real-valued response, quality is measured by mean squared error,

the average squared gap between what happens and what we predicted. We want it small.

We cannot compute it. The expectation is over the distribution of Y given X, which is the very thing we do not know. The obvious workaround is to approximate it with the data we have — the training set Tr = {(xᵢ, yᵢ)}ᵢ₌₁ⁿ:

And this is where it goes wrong. Minimising the training MSE on its own leads straight to overfitting: a function that passes exactly through every training point drives it to zero, a perfect score, while generalising terribly. The polynomial with too many powers and the neighbourhood shrunk too small both get there.

Return to the child learning what a tree is. How would you check that the learning actually happened, rather than five specific pictures being memorised? You would not show the same five pictures again. You would find a tree the child has never seen and ask.

That is the whole idea of a test set. Hold out data Te = {(xᵢ, yᵢ)}ᵢ₌ₙ₊₁ᵐ, never used in fitting, and evaluate on it:

Training error measures memorisation. Test error measures learning.

The behaviour of the two errors as flexibility increases is the central picture of the lecture. Take a fixed true f and fit three models of increasing flexibility — call them orange, blue, and green.

Read the three fits off those curves. Orange is inflexible and has high error on both training and test data: it underfits, and no amount of held-out data will rescue it because the model class simply cannot represent the truth. Green nearly interpolates, so its training error is tiny and its test error is large: it overfits. Blue does well on both and is the one to choose — and note that the training curve alone would never have told you that, since it keeps improving all the way to green.

How much flexibility is right depends on something we cannot see: the true f. When f is smooth, an inflexible model close to a straight line already does well, and adding flexibility buys almost nothing.

When f is less smooth — wiggly, with structure at a finer scale — the story reverses. The inflexible linear model performs poorly, and more flexible fits cut test error substantially at comparatively small cost.

So the right amount of flexibility depends on how complex the true f actually is, which is never known in advance. That is not a defeat — it is the reason model selection exists as a discipline, and the reason test data are not optional.

What the error is made of: bias, variance, and the irreducible part

The U-shaped test curve is an empirical fact. The next question is why it is shaped that way, and the answer comes from decomposing the error exactly.

Fix the setup. We have f̂ trained on training data, and a new point (x*, y*) with y* = f(x*) + ε*. The trick is to add and subtract 𝔼[f̂(x*)] — the average prediction our method would make at x* across different training sets — inside the square. Everything falls out into three pieces:

where

A word on what the expectation runs over, because it is easy to get this wrong. It is over y* and over the randomness in the training responses yᵢ — that is, over ε₁, …, εₙ, ε*. The feature values xᵢ and x* are held fixed throughout this course.

Now read the three terms.

𝕍[ε] is the irreducible error. It does not depend on f̂ in any way. Choose a better method, gather more data, think harder — you still pay it. It reflects the intrinsic difficulty of the problem, the part of Y that X simply does not determine.

The other two terms form the reducible error, and these are the ones our choices affect.

Variance, 𝕍[f̂(x*)], measures how much the prediction would change if we had drawn a different training set. A method that sticks closely to its training data has high variance: shift the data slightly and the fitted function moves a lot.

Bias, bias(f̂(x*)), measures the error from approximating f with an f̂ that structurally cannot capture its nuances. A linear method estimating a genuinely non-linear f has high bias — no training set will fix it, because the plane cannot bend.

Flexibility does not reduce error. It moves error from bias into variance.

That is the mechanism behind the U-shape. More flexible methods have lower bias and higher variance; less flexible methods have higher bias and lower variance. The test curve falls while bias reduction dominates and rises once variance takes over. And it is why no single method is superior to every other — the right balance depends on the true f, which nobody gets to see.

Classification: the Bayes classifier and k-nearest neighbours

Everything so far assumed Y ∈ ℝ. Plenty of problems are not like that. Tomorrow’s weather is rain, cloudy, or sun. A tumour is one of several cancer types. A recommender picks one of 328 films. These are qualitative responses, and their categories generally have no natural ordering — there is no sense in which cloudy is halfway between rain and sun. The lecture works in the binary case, Y ∈ {yes, no}, which we identify with Y ∈ {0, 1}.

Mean squared error stops being meaningful here, and the reason is worth stating plainly: the coding is arbitrary. Label the classes {0, 1} and you get one error magnitude; label them {0, 10} and you get a different one, for the identical set of mistakes. A measure that depends on an arbitrary choice of labels is not measuring the problem.

What does make sense is counting mistakes. The error rate is

the fraction of observations classified incorrectly. The same train/test split applies as before: the test error rate

approximates the population quantity

the probability of misclassifying a new point.

Now we can ask the question that produced the conditional expectation in the regression case: what is the best possible rule? Conditioning on X = x,

so minimising the error rate is exactly the same as maximising the probability of being right. If we knew the distribution of Y given X, the ideal rule would be

that is: predict 1 when P[Y = 1 | X = x] > 0.5, and 0 otherwise. This is the Bayes classifier, and it occupies the same position in classification that the conditional expectation occupies in regression — theoretically optimal, generalising cleanly to more than two classes, and completely inaccessible, because we do not know Y | X.

Its error rate is

which is the irreducible error of the classification problem. Even the best possible classifier misclassifies at this rate, because the label carries intrinsic randomness given X.

Inaccessible targets call for estimates, and the estimate here is the classification cousin of local averaging. Replace the unknown P[Y = j | X = x] with the proportion of nearby training points in class j:

where N holds the k training points closest to x. Which is to say: take a majority vote among the k nearest neighbours.

It works remarkably well. With k = 10 on the lecture’s simulated example, the kNN boundary closely resembles the Bayes boundary — a rule built from nothing but local vote counts recovering the optimal rule’s shape.

And the choice of k reproduces exactly the trade-off we have been circling all along. With k = 1, every region of the plane is decided by a single nearest point: the boundary becomes wildly irregular and unstable, which is high variance and overfitting. With k = 100, the vote is taken over so wide an area that the boundary smooths out and misses the true curvature — high bias, underfitting.

The parameter has a different name and the loss has a different form, but the tension is the one from the bias-variance decomposition, unchanged.

Conclusion

The lecture covers a lot of ground and it is really one continuous argument.

Machine learning inverts classical programming: rather than writing rules that turn data into answers, we supply data and answers and learn the rules. Made precise, the rule we are after is the conditional expectation f(x) = 𝔼[Y | X = x] — best predictor under squared error, whether justified by an additive-noise model or by the choice of loss alone — and it is unknowable, so it must be estimated. We estimate it either by assuming a form and fitting parameters, or by declining to assume and letting local data decide, and either strategy can be pushed until it merely memorises the training set. Which is why accuracy is judged on held-out data, where the test error curve is U-shaped, and why that shape has an exact explanation: expected error at a point is irreducible noise plus squared bias plus variance, with flexibility trading the second against the third. When the response is a category rather than a number, the loss becomes the error rate, the optimal rule becomes the Bayes classifier, kNN estimates it from data, and the choice of k recreates the same trade-off.

Two ideas underlie the whole course from here. The thing we want is a conditional feature of the distribution of Y given X — an expectation for regression, a maximiser for classification — and it is never directly available. And every method for approximating it sits somewhere on a spectrum trading structure against flexibility, bias against variance, prediction against interpretation.

The lectures ahead spend their first half in the interpretable, linear-based world, extending it with splines and generalized linear models, before moving to the flexible and far more opaque methods that define the current state of the art.

Thanks for reading! If you found this useful, consider following for more articles on machine learning and statistical modelling.

Using ten years of Rio de Janeiro temperatures to hunt for daily, weekly, and annual patterns — one chart at a time.

Introduction

There’s a temptation in time series analysis to jump straight to decomposition, modelling, and forecasting. I’ve done it. The problem is that if you skip the exploratory work — or rush it with a single time-series plot at one resolution — you can easily miss entire seasonal structures hiding in plain sight. As my first time series professor, Cris, once put it: for seasonal series, get the seasonal structure right and even the most naive model can surprise you.

This article is a visual walkthrough of that exploration. We use ten years of 2 m air temperature in Rio de Janeiro 2015–2024, Open-Meteo ERA5 — [Historical Weather API docs) as a running example. Rio is instructive: the physics are clean, the signal is strong, and we have the same variable at three granularities — hourly, daily, and monthly — which means we can go looking for daily, weekly, and annual seasonality with the right tools for each.

The goal is not a gallery of charts. Each figure is a different way of asking the same few questions: does temperature repeat within the day (across hours), within the week (across weekdays), and within the year (across months and seasons)? We stay visual throughout; formal tests and models come later.

Multiple Seasonalities: A Quick Framing

First things first, what is seasonality? Let me give you two answers — one intuitive, one more precise.

Intuitively, seasonality is just a pattern that repeats on a predictable calendar rhythm: hotter summers, quieter Sundays, busier Decembers. More precisely, a time series is seasonal when it exhibits periodic fluctuations with a fixed and known period m — the number of observations in one full repetition of the pattern. Both definitions point to the same idea: structure you can anticipate, if you know where to look. It is important to say that periodic pattern does not mean a deterministic pattern. The seasonal component itself can, and commonly is, stochastic — aka random — but the period m must be fixed.

In Forecasting: Principles and Practice (3rd ed.), Hyndman and Athanasopoulos describe multiple seasonal patterns as the presence of more than one repeating seasonal structure in the same series — and they note it is routine in high-frequency data. Their electricity-demand example from Australia shows all three at once: lower use overnight, different weekday vs. weekend behaviour, and higher demand in summer and winter. Surface temperature in Rio is a different variable, but the logic is the same. We should look for daily, weekly, and annual seasonality, and we should expect some to be strong and others simply absent.

One naming point matters before we go further. In FPP’s sense, daily seasonality does not mean “the whole series wiggles every calendar day in the long-run plot.” It means the pattern across hours within a day — cool before sunrise, warm in early afternoon — that repeats day after day. That is the profile you can only see clearly when you work with hourly data or hour-of-day summaries.

The seasonal period m counts how many observations fit in one full repeat of a pattern. Crucially, m depends on the granularity of the series — we are not changing the physics when we switch between CSVs, only changing how many data points fall inside one seasonal period:

The dataset

Three series from the same ERA5 reanalysis, all measuring `temperature_2m_c`:

We move from wide views of the full series toward plots that isolate each seasonal scale. Heatmaps show means across two calendar dimensions; boxplots add spread by hour, weekday, and month; profile plots show the typical seasonal shape and how much it varies from year to year or day to day.

Setup and Reproducibility

All the code and data download scripts live in the companion repository on GitHub: SeasonalityGraphs project folder. The README walks you through the environment setup, how to fetch the Open-Meteo ERA5 data, and how to regenerate all 18 figures in one command.

First Look: Three Resolutions

The natural starting point is to plot the entire series at three granularities — hourly on top, daily in the middle, monthly at the bottom.

Annual seasonality jumps out immediately. On the daily series (m = 365), you can count ten warm–cool swings — austral summers (December–February) high, June–August low. The monthly series (m = 12) strips almost everything else away and leaves that annual rhythm alone. You can see that even though the repetition is clear, the seasonal pattern is not deterministic; it is stochastic. This means that every summer is warmer then every winter, but how much warmer changes over the years.

Daily seasonality (m = 24, hourly data only) does not show up here. The top panel looks like textured noise at this zoom; you cannot read the hour-of-day profile from a decade-long hourly trace. That is the main lesson of Fig. 01: a single time-series plot at full length is well-tuned for annual structure, and poorly suited for everything else. To investigate daily seasonality you need plots that explicitly aggregate by hour — heatmaps, boxplots by hour, and hourly profiles. Let’s build those next.

Two Calendar Dimensions at Once

A line plot shows when temperature moves. A heatmap shows two calendar axes simultaneously: the mean temperature in each cell. On hourly data, that is where annual and daily seasonality appear side by side — and where you start to notice that their interaction is not trivial.

Hour × Month — Fig. 02

Read along a horizontal band at a fixed hour: temperature rises and falls with the month — annual seasonality on hourly data (here collapsed to month-of-year for readability). Read up a column at a fixed month: cool before sunrise, peak around noon–2 pm, cooler in the evening — daily seasonality (m = 24).

Now compare January and July. In January, the gap between the 6 am row and the noon row is large; in July, it is noticeably smaller. The shape of the daily cycle is not the same in every month — hour and month interact in a way that a simple additive decomposition will not capture perfectly. Worth keeping in mind when you move to modelling.

Hour × Weekday — Fig. 03

We include weekly seasonality (m = 168 on hourly data) for completeness, but the answer is essentially foregone: there is no physical reason Sunday should be systematically warmer than Wednesday. The atmosphere does not follow the calendar week. The heatmap matches that expectation — all seven columns show the same gradient from cool pre-dawn to warm midday. Nothing shifts from Monday to Sunday. Weekly seasonality is absent, and we will see the same flat story repeated in every subsequent plot that touches weekdays.

Seasonal Levels and Spread by Calendar Feature

Heatmaps show means only. A warm average at noon could mean every day is consistently warm at noon, or it could hide a mix of very hot and mild days. Boxplots let us read both the typical level (median) and the spread (IQR, whiskers) by calendar category — the same three seasonal questions, now with an extra layer of information about variability.

By Hour of Day — Fig. 04

Daily seasonality is strong and easy to read. Early-morning hours (2–6 am) have tight boxes — consistently cool. Midday hours (11–15) have wider IQRs and longer upper whiskers — not only warmer on average, but more variable from day to day. The spread is itself a feature of the daily cycle, not just the level. Remenber that those boxplots are showing the hourly distribution of the full years. One interesting extra step could be to divide this plot into 4 seasons to answer the question: does the daily variability changes across seasons? Perhaps during the summer the days are consistenly warmer. Perhaps winter temperatures changes a lot during the day.

By Weekday — Fig. 05

The expected null: no weekly seasonality. All seven days sit on top of each other. A weekday dummy in a temperature model would add nothing.

By Calendar Month — Fig. 06

The median steps smoothly through the year — classic southern-hemisphere annual pattern. Summer months have wider IQRs than mid-winter; Rio is not only hotter in summer but also more variable day to day. Winter months (especially July–August) are narrower and more predictable. That asymmetry in spread matters: it shows up later in the profile plots too.

Typical Seasonal Profiles

Now we plot the average seasonal shape for each scale and overlay many realisations — one line per year, or per day — along with a P10–P90 envelope and a red mean line.

A wide envelope means the seasonal pattern is stable in shape but not in level from one year or day to the next.

Annual Profile by Month — Fig. 07

Built from hourly data aggregated to mean temperature by calendar month. This is a clean picture of annual seasonality — m = 12 at monthly resolution. The ten years are nearly parallel: same timing, similar amplitude. The summer peak is where years diverge most, consistent with the boxplot in Fig. 06. Notice that the spread of the P10-P90 interval follows perfectly the spread of the boxplots in the previous figure (look at March and October for a clear comparison).

Daily, Weekly, and Annual Profiles — Fig. 08

Top panel — daily seasonality (m = 24): Mean temperature rises from ~22°C pre-dawn to ~26–27°C between 11 am and 2 pm, then falls. The envelope is tight at 6 am and fans out by noon (nearly 15°C between P10 and P90) — many days share the same timing of the peak even when they disagree on how hot the afternoon gets.

Middle panel — weekly seasonality: The mean is a flat line. The envelope is two horizontal bands. No day of the week stands out. This panel is the formal version of what we already called obvious.

Bottom panel — annual seasonality (m = 365 on daily data): Ten years on the day-of-year axis follow the same broad arc — hot at the year edges (austral summer), cool around days 190–210 (winter) — but they cross more than the monthly profile in Fig. 07 does, because daily values carry more within-year weather noise. Year-to-year differences become hard to separate from day-to-day variability; still, every year follows the same gross shape. Some years run hotter or cooler in particular summers or winters, but the seasonal clock stays put. If you wish to explore long term patterns, such as temperature changes over the year, maybe a high frequency data is not ideal, since it may introduce “noise” from other high frequency patterns.

Does the Annual Pattern Repeat Across Years?

The bottom panel of Fig. 08 already shows the annual structure on daily data. Here we zoom in with a cleaner view — one point per month, one line per year — to get a sharper picture of year-to-year consistency.

Monthly Resolution — Fig. 09

With monthly aggregation (m = 12), the picture is clean. All ten years follow the same annual rhythm. Winters cluster closely; Jan–Feb separate more across years (~4°C between the warmest and coolest in this window). The divergence is mainly in amplitude at the summer peak — not in when the trough arrives. July or early August stays the coldest slice across all years.

Zooming Inside Calendar Strata

Subseries plots fix one calendar bucket — a month, or a weekday — and show the series through time inside that stratum. They answer a natural follow-up question: once we know the global annual or weekly structure, what does a single slice actually look like up close?

Twelve Months — Fig. 10

Each panel is a slice of annual seasonality: January’s mean sits well above July’s. Within a given month, points scatter around the dashed mean without a clear decade-long slope in most panels — but January and other summer months span a noticeable range between years. In this daily-based view, that spread is already clear: the summer months are visibly wider year to year than winter, and January varies by roughly 4°C across the 2015–2024 window. Modelling “January in Rio” as a fixed seasonal level plus noise would miss real interannual movement in the warm season.

Seven Weekdays — Fig. 11

The clearest restatement of the obvious: Monday, Wednesday, and Saturday share the same mean level and the same decade of annual ups and downs. Weekly seasonality has nothing to add here. Calendar weekday dummies would be dead weight for this variable.

What We Learned

Rio’s 2 m temperature from 2015–2024 is a clean example of multiple seasonalities with very unequal strength.

Annual seasonality dominates the narrative whenever we step back from hourly noise — roughly 5–6°C from the winter median to the summer median, stable in timing across years but with meaningful spread in hot months (up to ~4°C between years in Jan–Feb). On daily data m = 365; on monthly data m = 12.

Daily seasonality is equally real but invisible on a long hourly plot (Fig. 01). It appears once we aggregate by hour: a ~4–5°C swing from pre-dawn to afternoon in the mean profile, with a wider envelope at midday and a larger daily amplitude in summer than in winter (Figs. 02, 04, 08). On hourly data m = 24.

Weekly seasonality is absent — as physics would predict. Heatmaps, boxplots, profiles, and subseries all tell the same flat story. A m = 168 or m = 7 term would not earn a place in a model here.

Spread matters too: afternoons and summer months are more variable, not just warmer. Hour and month interact slightly in the daily shape (Fig. 02). None of this replaces formal inference, but it narrows what a forecaster should plan for: at least daily and annual seasonal structure on hourly data, no weekly term, and interannual variation in summer levels that is worth modelling explicitly.

Investigating multiple seasonalities takes multiple visualization techniques. A plain time-series plot of the full series is the one view most likely to hide high-frequency structure — and the visual noise in dense hourly data can make us overlook daily or weekly patterns altogether.

Next steps

This article stays exploratory and graphical. Natural follow-ups — aligned with FPP3 Chapter 12 on complex seasonality — include:

• MSTL / multiple-seasonal decomposition on hourly data to separate daily and annual components simultaneously
• ACF at seasonal lags — 24, 168, and 8766 for hourly series; 7 and 365 for daily
• Dynamic harmonic regression (Fourier terms) or TBATS for multiple seasonal periods

The next piece in this series moves from “we can see the structure” to “we can estimate and forecast it.”

• Code and data:

MathNogMediumArticles/TimeSeriesEDA/SeasonalityGraphs at main · MathNog/MathNogMediumArticles

• Data: Open-Meteo ERA5

🏛️ Historical Weather API | Open-Meteo.com

• Forecasting Principles and Practice

Forecasting: Principles and Practice (3rd ed)

Thanks for reading! If you found this useful, consider following for more articles on time series analysis and forecasting.

Before learning which model to use, it helps to understand what a model actually is.

The World Is Complicated

The world is complicated — and we can try to interpret reality in many ways.

One approach is to take a slice of reality and represent it through a set of assumptions expressed in the language of mathematics. Think about it: if we want to understand how the average daily temperature of a city evolves, we can suppose that today’s temperature shouldn’t differ dramatically from yesterday’s. Nearby days should carry information about the next day. But we also know that average temperatures shift with the seasons, so to predict tomorrow’s temperature well, it’s also informative to look at what temperatures looked like at this time last year. That means we might need to reach back at least 12 months into the past.

This kind of reasoning — laying down assumptions, translating them into equations, then checking if they were reasonable — is exactly what modelling is.

This article is about the foundations. Not which model you should use, but what a model is, why the random error term matters more than people give it credit for, and how to tell if the model you chose was a good idea. By the end, you’ll be able to look at a model equation and understand what every piece of it is doing.

Let’s dive in 🚀

What Is a Time Series?

Let me give you two answers — one intuitive, one slightly more precise.

Intuitively, a time series is a numerical measurement (temperature, price, number of passengers, electricity demand…) recorded at regular time intervals (every minute, every hour, every month…) that tends to exhibit some degree of randomness.

More precisely, a time series is a realisation of a discrete-time stochastic process with a continuous state space. If that sounds like a mouthful, don’t worry — the intuitive definition will carry us through this article. The formal version just says that the sequence of values you observe is one possible outcome drawn from a random process that unfolds in time.

One of the most iconic examples is the AirPassengers series: monthly counts of international airline passengers from 1949 to 1960. It shows up in practically every introductory time series course, and for good reason — it is simple, clean, and has two features that make it interesting: an upward trend and a repeating seasonal pattern.

import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
df = pd.read_csv("data/airline.csv", parse_dates=["date"])
df = df.sort_values("date").reset_index(drop=True)
fig, ax = plt.subplots(figsize=(9, 4))
ax.plot(df["date"], df["passengers"], color="#2B2B2B", linewidth=1.6)
ax.set_title("Monthly International Airline Passengers (1949–1960)")
ax.set_ylabel("Passengers (thousands)")
ax.xaxis.set_major_formatter(mdates.DateFormatter("%Y"))
plt.tight_layout()
plt.show()

We’ll come back to this series throughout the article.

What Does It Mean to Model a Time Series?

Modelling a time series means finding a function that, given the past values of the series (and possibly other inputs), produces a useful approximation of the current or future value. Our goal is to describe, through one or more equations, how the series behaves given its own history:

Here, yₜ is the value at time t, f is the functional form we choose, and θ is the vector of unknown parameters the model needs to be estimated.

But — and this matters — it would be arrogant to assume any model could perfectly capture reality. So we add a random error term:

This εₜ is not just a cosmetic addition. It is the part of reality the model will not (and should not be expected to) explain. It is what separates a forecasting method from a statistical model.

Once we write down εₜ, we have to make assumptions about it. The most important: E[εₜ] = 0 (zero mean). This says the model shouldn’t be systematically wrong — it shouldn’t consistently overestimate or underestimate. Other common assumptions include no autocorrelation (the errors at different time points are uncorrelated) and homoscedasticity (constant variance over time). For now, zero mean is the one we’ll keep an eye on.

A model without an explicitly defined error term is just a forecasting heuristic. The error term is what makes it a statistical model.

Because εₜ is, by assumption, random and structureless, it should look like white noise — no patterns, no trend, no seasonality. This becomes the key diagnostic tool we’ll use at the end.

A Ladder of Models

Now that we have the general equation, let’s climb a ladder of increasingly rich functional forms — each one a different assumption about how the past informs the future.

The Random Walk

The simplest possible model says that the best guess for today is just yesterday’s value:

This is the Random Walk — arguably the most famous model in time series analysis. It shows up in finance, physics, and biology. It makes no attempt to explain why the series moves; it just says movement is random around wherever we were.

A Moving Average of the Past

One step up: what if the best guess is the average of the last three observations?

This is more stable than the Random Walk (outliers get averaged out), but it still has no parameters to estimate — the weights are fixed at 1/3 each.

The AR(p) Model

Now let’s let the data decide how much each past value matters. Instead of fixed equal weights, we estimate a coefficient for each lag:

This is the AR(p) model — Autoregressive of order p — and it is one of the most important models in classical time series. The β coefficients are unknown and must be estimated from the data. In machine learning language, that’s what people call “training the model.”

The two simulations below illustrate how the choice of functional form shapes the behaviour of the series — even when both are driven by identical noise.

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)
n_total = 500   # simulate long enough to burn in initial conditions
n_plot = 200    # plot only the last segment
phi = 0.85

eps = np.random.normal(0, 1, n_total)
rw = np.zeros(n_total)
ar1 = np.zeros(n_total)
rw[0] = eps[0]
ar1[0] = eps[0]

for t in range(1, n_total):
    rw[t] = rw[t - 1] + eps[t]              # y_t = y_{t-1} + eps_t
    ar1[t] = phi * ar1[t - 1] + eps[t]      # y_t = 0.85 * y_{t-1} + eps_t

rw = rw[-n_plot:]
ar1 = ar1[-n_plot:]

fig, axes = plt.subplots(1, 2, figsize=(11, 4))
axes[0].plot(rw, color="#2B2B2B", linewidth=1.4)
axes[0].set_title("Random Walk\n$y_t = y_{t-1} + \\varepsilon_t$")
axes[1].plot(ar1, color="#2B7BB9", linewidth=1.4)
axes[1].set_title("AR(1)\n$y_t = 0.85\\,y_{t-1} + \\varepsilon_t$")
for ax in axes:
    ax.set_xlabel("Time")
plt.tight_layout()
plt.show()

This is the key insight behind the ladder metaphor: each model is just a different choice of f. The more structure we build in, the more assumptions we’re making — and the more we can be wrong.

The landscape of available models is vast. To give you a sense of the space:

- Classical statistical models — SARIMA family (MA, AR, ARMA, ARIMA, SARIMA), Unobserved Components models, Exponential Smoothing (ETS)

- Machine learning models — Linear/Ridge/Lasso regression, SVR, Random Forests, Gradient Boosting (XGBoost, LightGBM, CatBoost)

- Deep learning models — RNN, LSTM, GRU, CNN, Transformers

How Do You Choose a Model?

That choice depends mainly on two things: how the series behaves, and what the problem demands.

The behaviour of the series is uncovered through exploratory data analysis. Plot the series, look for trend and seasonality, check whether it looks stationary, inspect autocorrelation — the usual toolkit. We’ll dig deeper into that side of things in future articles in this series.

The problem context shapes how complex and interpretable the model should be. If you need fast, repeated re-estimation on a tight compute budget, simpler classical models are often the right tool. The same goes when you need to explain why a forecast looks the way it does. If the application only cares about forecast accuracy and you have enough data to train it, a non-interpretable black-box from the ML or deep learning shelf can be a strong option — as long as you’re comfortable trading transparency for flexibility.

There is no single winner in the list above — only models that fit the data and the job better or worse.

For the rest of this article, we’ll stay in the classical territory and look at one family: ETS. If you wish to dive into time series modelling and forecasting the best approach is to follow Rob Hyndman free book Forecasting Principles and Practice.

The ETS Family

The ETS family — standing for Error, Trend, Seasonality — grew out of exponential smoothing methods developed in the late 1950s (Brown, Holt, Winters). As described in Forecasting: Principles and Practice, the Pythonic Way Chapter 8, exponential smoothing produces forecasts as weighted averages of past observations, where the weights decay exponentially — more recent observations count more.

What makes the ETS family powerful is its ability to combine three components in different ways:

- Error (E): additive (A) or multiplicative (M)

- Trend (T): none (N), additive (A), or multiplicative (M)

- Seasonality (S): none (N), additive (A), or multiplicative (M)

Each combination gives a different model. The taxonomy below shows all possible combinations:

A few members of this family have famous names:

- ETS(A,N,N) — Simple Exponential Smoothing: no trend, no seasonality

- ETS(A,A,N) — Holt’s linear method: additive trend, no seasonality

- ETS(A,A,A) — Holt-Winters additive: additive trend and additive seasonality

Look back at Figure 1. The airline series has an obvious upward trend and clear seasonal swings every 12 months. Intuitively, the ETS(A,A,A) should handle that well. The ETS(A,N,N) — which knows nothing about trends or seasons — should struggle.

Let’s see this in practice. We fit both models with statsmodels’ ETSModel, which implements the state-space ETS framework (see the official ETS examples) and estimates smoothing parameters plus initial states by maximum likelihood:

import pandas as pd
from statsmodels.tsa.exponential_smoothing.ets import ETSModel

df = pd.read_csv("data/airline.csv", parse_dates=["date"])
y = df["passengers"]

# ETS(A,N,N): additive errors, no trend, no seasonality
fit_ann = ETSModel(
    y, error="add", trend=None, seasonal=None
).fit(maxiter=10000, disp=False)

# ETS(A,A,A): additive trend and additive seasonality (period 12)
fit_aaa = ETSModel(
    y,
    error="add",
    trend="add",
    seasonal="add",
    seasonal_periods=12,
).fit(maxiter=10000, disp=False)

fitted_ann = fit_ann.fittedvalues
res_ann = fit_ann.resid
fitted_aaa = fit_aaa.fittedvalues
res_aaa = fit_aaa.resid

# Smoothing parameters (α, β, γ) — see fit.summary() for full output
params_ann = dict(zip(fit_ann.model.param_names, fit_ann.params))
params_aaa = dict(zip(fit_aaa.model.param_names, fit_aaa.params))

The companion script `code/what_model_means_models.py` wraps the same calls in `ets_ann_fit()` and `ets_aaa_fit()` for the figures and CSV exports.

After fitting, we get:

The residual standard deviation is already telling us something: ETS(A,A,A) leaves much less unexplained noise than ETS(A,N,N). The fitted values make it even clearer. The fitted values are the outputs of the fitted model, given its inputs — aka the past values of the series and the estimated parameters.

The visual speaks for itself — the well-specified model hugs the data; the misspecified one trails behind, forever chasing a trend it never learned.

How Do We Know If the Model Is Good? Residual Diagnostics

Now we return to εₜ. Once a model is estimated, we can compute the residuals:

These are the realisations of εₜ under our estimated parameters — what the model got wrong at each time step. And here’s the key idea: if the model was well specified, the residuals should look like random noise.

Remember our assumptions about εₜ? They were:

1. Zero mean — the model is unbiased

2. No temporal structure — no patterns left unexplained

We can check both visually by simply plotting the residuals over time. If we see patterns — a trend, seasonal waves, oscillations — we know the model left structure on the table. Below, `res_ann` and `res_aaa` come from `fit_ann.resid` and `fit_aaa.resid` after the statsmodels fits above.

fig, axes = plt.subplots(2, 1, figsize=(10, 6), sharex=True)

for ax, residuals, label, color in zip(
    axes,
    [res_ann, res_aaa],
    ["ETS(A,N,N) residuals", "ETS(A,A,A) residuals"],
    ["#E05C2B", "#2B7BB9"],
):
    ax.plot(df["date"], residuals, color=color, linewidth=1.3)
    ax.axhline(0, color="#2B2B2B", linewidth=0.9, linestyle="--")
    mean_val = residuals.mean()
    ax.axhline(mean_val, color=color, linewidth=1.1, linestyle=":",
               label=f"mean = {mean_val:.1f}")
    ax.set_title(label)
    ax.legend()

plt.tight_layout()
plt.show()

The ETS(A,N,N) residuals are practically screaming seasonality. Peaks and troughs repeat year after year — clear evidence the model missed something important. The ETS(A,A,A) residuals, on the other hand, look much more like random noise. Their mean is close to zero (≈0.0 vs. ≈2.2 for the misspecified model), and no obvious temporal pattern emerges.

This is residual diagnostics in its simplest form. More formal tools exist — autocorrelation functions (ACF/PACF), Ljung-Box tests, heteroscedasticity checks — but even a visual inspection of the residual time plot can tell you a great deal.

Putting It All Together

Let me summarise what we’ve built:

Modelling a time series means choosing the functional form of f — which lags to include, how to combine them, whether to account for trend and seasonality — and then estimating any unknown parameters θ from the data (here via statsmodels’ likelihood-based `ETSModel.fit`). Once estimated, the model’s assumptions live in the residuals: if εₜ was assumed to be structureless, the residuals should look structureless too.

The Random Walk is the simplest possible f. The AR(p) is a flexible linear combination of past values. The ETS family goes further, decomposing the series into error, trend, and seasonal components that interact — either additively or multiplicatively.

Choosing the right model depends on two things: the behaviour of the series (visible in exploratory data analysis) and the context of the problem (how much interpretability you need, how fast re-estimation must be, how much data is available for training). There is no universally best model — only models that are more or less compatible with the data and the task at hand.

Future articles in this series will cover the specific models, time series EDA and eventually the machine learning approaches that make fewer assumptions about the functional form of f. Stay tuned 🥁

All code in this article is available in my GitHub account. Feel free to run everything locally — install dependencies from `requirements.txt` (NumPy, pandas, matplotlib, statsmodels).

If you enjoyed this article, follow for more on time series, forecasting, and mathematical modelling. And if something here helped you understand something that was once murky — that’s exactly why I write these.

In this short article, you will learn the basics of a super flexible and general time series model framework called Score Driven Models. They were developed simultaneously by Harvey [2013] and Creal et al. [2013]. Known by different names — generalized autoregressive score (GAS), dynamic conditional score models, or simply score-driven models, they offer several advantages for time series modeling, as described below. Firstly, these models can handle non-Gaussian time series, where the conditional probability distribution is not normal. This capability is crucial for cases where the normal distribution is inadequate, such as wind speed series (which only take positive values) and financial return series (which exhibit excess kurtosis). Count data, such as the number of items sold in a store, also benefit from this flexibility, as they require discrete distributions. The ability to model non-normal data makes GAS models highly applicable across various domains.

In this article, we will use a Natural Inflow Energy time series from Brazil, which is naturally non-gaussian, to exemplify the use of this type of model with an open-source package written in Julia called ScoreDrivenModels.jl.

In order to keep this article short, I will not dive into the full length of the theory behind score driven models. To do so, I strongly recommend diving into the following resources:

• Gas Model Website with a lot of published papers and material regarding score-driven models;
• Harvey original paper with the original description of the Dynamic Conditional Score model;
• Cral et al. original paper with the original implementation of the GAS(p, q) model.
• My undergraduate final project, submitted as a complete paper at the SBPO (Brazilian Symposium on Operations Research)

What is a Score Driven Model?

The best way to introduce a score-driven model is to define the fundamental steps that we need to follow to model a time series with this framework. Since they are a probabilistic model, out first step is to define a conditional probability distribution to be followed by the time series. Note that this step is super general! We can choose any distribution that we think will best describe the nature of the time series. Gaussian, Lognormal, Student-t, F, Gamma, and so on.

This definition is saying the out time series y follows a conditional probability distribution of both the past observations of the series, static parameters and time-varying parameters (which are of crucial importance to this framework). One this distribution is defined, we need to choose how the time-varying parameters will evolve over time, since they are (who would guess?) time-varying. There are two main options: (a) an ARMA(p, q) process, following Creal et al., giving origin to a GAS(p,q) model; or (b) an unobserved components framework, following Harvey’s proposition. We will choose the first one, since it is the one implemented inside the ScoreDrivenModels.jl package.

The GAS(p,q) model says that the time-varying parameters follow an ARMA(p,q) process. Yes! the same from the ARIMA family! Actually, the GAS(p,q) model generalizes both ARMA and GARCH models. Since we are imposing an ARMA(p, q) process in the time-varying parameters, this makes the GAS(p,q) framework inappropriate for non-stationary time series. This is a core point of this methodology! The unobserved components framework, on the other hand, is appropriate for non-stationary time series! Let’s define an ARMA(p, q) process to the time-varying parameters dynamic.

Wait a minute! There is a new symbol in this definition which we have not presented yet: s. This is, indeed, the core of the proposed framework. It is the standardized score of the time series. Let’s define, first, who is the score: it is the first derivative of the log-likelihood function with respect to the parameter of interest. This has some intuition: is says how much improvement we have on the likelihood of our model if we change the parameter of interest (mean, variance, and so on). It is the score that drives the dymanic of the score-driven model!

Finally, we can define the standardized score, which is simply a scaled version of the score obtained by a (highly non trivial) multiplication. I will not dive into the properties of this scaled version, but it is a key element of the framework.

The last piece of theory I would like to talk about is the model’s estimation. As usual, it is done with maximum likelihood estimation, in which we aim to find the optimal parameters that maximize the likelihood of the model given the observed time series. This is particularly straightforward, provided an efficient nonlinear optimization algorithm is used, requiring only the implementation of the updating function and the evaluation of the likelihood function at a particular optimization point.

That’s it for theory! Let’s talk about the ScoreDrivenModels.jl package.

The ScoreDrivenModels.jl package

This package, developed by our laboratory LAMPS implements score driven models of type GAS(p,q) for a variety of available distributions with 1 or more time-varying parameters.

Take a look into the package documentation to a full report of it’s functionality! You can also dive into the package’s article.

The table below summarizes the available distributions where the identity scale means d=0, inverse scale means d=-1 and inverse square means d=0.5.

Moreover, the package also allows for new distributions to be implemented with pull requests. All one have to do is to define some values such as the score, Fisher information and link functions to the desired distribution.

There is a number of available optimizers to fit the model with MLE, such as Nelder Mead, LBFGS and IPNewton. However, there is a fundamental trick!

Those optimizers are unconstrained algorithms, which means that they do not guarantee any type of restrictions of the distribution parameters. This is a huge problem, since variance, for instance, can never be negative!

The trick to guarantee this type of restrictions is through link functions, which are described in the package documentation, here.

That’s enough for plain text and mathematical equations. Let’s dive into the action!

Modeling a non-gaussian time series with ScoreDrivenModels.jl

First of all, we need the following packages:

using Dates
using Plots
using DelimitedFiles
using Distributions
using ScoreDrivenModels

Loading a time series

This is a monthly time series of the natural inflow energy of the Northeast region of Brazil, available inside the package.

I have downloaded it from the package repo and saved it in the `data` folder in my GitHub repository. The following code loads the time series and divide into train and test sets.

dates = collect(Date(1961):Month(1):Date(2000, 12));
y = vec(readdlm("data/nie_northeastern.csv"));

H  = 60
T  = 240
last_obs  = length(y) - H
first_obs = last_obs - T

y_train     = y[first_obs:last_obs];
y_test      = y[last_obs+1:end];

dates_train = dates[first_obs:last_obs];
dates_test  = dates[last_obs+1:end];

We can also plot the time series to visually inspect it!

plot(dates_train, y_train, label = "Train")
plot!(dates_test, y_test, label = "Test") 
plot!(title="Brazil Northeastern Natural Inflow Energy", xlabel="Date", ylabel="Value")
plot!(xformatter = x -> Dates.format(Date(Dates.UTD(x)), "yyyy"))

As we can clearly see, the time series has no visual trend, mut a strong seasonal pattern. Moreover, this series does not follow a gaussian distribution, because it is impossible to have negative inflow!

Defining the model

We define the model using the ScoreDrivenModel function from the package. The function takes the following arguments:

- p_lags: The lags of the autoregressive part of the model.

- q_lags: The lags of the moving average part of the model.

- dist: The distribution of the model.

- d: The scale for the standardized socre of the model.

- time_varying_params: The parameters of the model that are time-varying.

This is a crucial part of the modeling process, as it is when, based on the time-series features, we define the model’s structure and the parameters that we want to estimate.

Since our time series is of monthly frequency and it is a time series of natural inflow energy, we expect to have a strong seasonal pattern (as shown in the previous plot). Because of that, we will add a seasonal lag of 12 months in both the autoregressive and moving average parts of the model.

Moreover, since natural inflow can never be negative, it is natural to use a distribution that is always positive, such as the log-normal distribution, from the Distributions.jl package. It is important to say that there is no technical reason not to use a Normal distribution. The model would be fitted and simulated without any issues, but we would possibly generate negative scenarios, which would not make sense in this context!

Finally, we will fit a Log-normal GAS model with both parameters as time varying, that is, both will follow an ARMA process with autoregressive and moving average lags of 1 and 12 months.

p_lags  = [1, 12]
q_lags  = [1, 12]
distrib = Distributions.LogNormal
d       = 0.0

gas = ScoreDrivenModel(p_lags, q_lags, distrib, d; time_varying_params = [1, 2]);

Estimating and Simulating the model

We can fit and forecast the model with two simple functions! Notice that we are doing a probabilistic forecast with 1000 scenarios!

fit_gas = fit!(gas, y_train);
forecast_gas = forecast(y_train, gas, test_months; S=1000);

Once the model is fitted and forecasted, we can visualize it’s results!

First, the fitted parameters.

We can also plot the fitted mean and both point and probabilistic forecasts.

point_forecast = forecast_gas.observation_forecast;
scenarios      = forecast_gas.observation_scenarios;
quantiles_obs  = forecast_gas.observation_quantiles;
fitted_values  = fitted_mean(gas, y_train);

plot(dates_train, y_train, label = "Train Values", color = :black)
plot!(dates_train, fitted_values, label = "Fitted Values", color = :blue)
plot!(dates_test, scenarios, label = "", color = :lightgray)
plot!(dates_test, y_test, label = "Test Values", color = :black)
plot!(dates_test, point_forecast, label = "Point Forecast", color = :red)
title!("Natural Inflow Energy Simulation from LogNormal GAS")
plot!(xformatter = x -> Dates.format(Date(Dates.UTD(x)), "yyyy"),
    legendcolumns = 4, legend = :outerbottom)

Even thought this may seen enough to present a time series modeling and forecasting package, there are some super interesting and powerful functionalities implemented that must be talked about!

Fitted and forecasted parameters

First, we can visualize each parameter fitted values and probabilistic forecast! Yes, even thought variance is, essentially, unobserved, we can visualize it with this package as follows.

params_fitted = score_driven_recursion(gas, y_train)
params_point_forecast = forecast_gas.parameter_forecast; #H x P
params_scenarios      = forecast_gas.parameter_scenarios; #H x P x S

for i in 1:size(params_fitted, 2)
 plot(dates_train, params_fitted[2:end, i], label = "Fitted Parameters", color = :blue)
 plot!(dates_test, params_scenarios[:, i, :], label = "", color = :lightgray)
 plot!(dates_test, params_point_forecast[:, i], label = "Point Forecast", color = :red)
    title!("Parameter $i for LogNormal GAS")
end

Residuals! Always check the residuals!

One crucial step in any time-series model workflow is to check the residuals of the model. Features like zero mean, homoscedasticity and lack or auto-correlation are essential to validate the fitted model and to allow ourselves to look into (and trust!) the out of sample forecasts. Even though the package provides a residuals diagnostics plot, it is not perfect, and I highly encourage you to implement your own diagnostics. We can access the simple one provided by the following code.

plot(fit_gas, size=(1000, 1000), title = "Residuals Diagnostics")

Which generates this result.

We can also compute both Pearson and Quantile Residuals automatically with the package so that we are able to, at least, proceed with a visual inspection. Allow me to simply paste those residuals definitions here, without a lot of explanation…

q_residuals   = quantile_residuals(y_train, gas);
std_residuals = pearson_residuals(y_train, gas);

Built-in Time Series Cross Validation

There is a built-in function for time series cross validation, in which we simply define the forecast horizon and the first timestamp index to fit the model. This function automatically computes the number of windows to fit and forecast the model, while returning different accuracy metrics across all windows and horizons. This is an excellent feature to access the proposed model consistency thought different sets of training and testing, improving our confidence in the model results.

first_index = 100
gas_cv = cross_validation(gas, y_train, H, first_index);

cv_abs_errors = gas_cv.abs_errors
cv_mean_crps  = gas_cv.mean_crps
cv_mae        =  gas_cv.mae

H, W = size(cv_abs_errors)

plot(1:H, cv_mean_crps, dims=2, label="Mean CRPS")
plot!(1:H, cv_mae, label="MAE")
plot!(xlabel = "Window", ylabel = "Error", title = "Mean Absolute Error, Mean CRPS and MAE across all $W windows by $H steps ahead",
titlefontsize = 11)

Thanks for reading this far! I hope this tutorial helped you to understand the core power and flexibility of ScoreDrivenModels to model and forecast non-gaussian (and gaussian, why not?) time series in Julia. Feel free to use the package yourself and let us know about your experience!

All the code used here is 100% written by me and available on my GitHub. The text, on the other hand, was written with the help of AI — because good tooling helps with writing too 😉

If you liked this package, be sure to check out the other time series tools we’ve built at LAMPSPUC. And stay tuned — more tutorials are on the way!

GitHub - MathNog/MathNogMediumArticles: Public Repository for my Medium

Who am I?

Hey there! I’m Matheus Nogueira, a Brazilian master’s student in Operations Research at LAMPS, an Optimization and Statistics Lab at PUC-Rio. I did my undergrad in Computer Engineering, also at PUC-Rio, in the Informatics’ Department.

Most of my work for the last years is about making sense of time series — mainly trying to predict the future with forecasting and simulation and I would love to share some interesting tools and concepts with you.

What’s LAMPS?

LAMPS is a research lab that mixes optimization, statistics, and machine learning — with people from both Electrical and Industrial Engineering at PUC-Rio. Besides academic research, we also work with companies and develop open-source packages, mostly in Julia.

If you’re into time series, you’ll probably like our GitHub. We’ve got packages for ARIMA, state space models, score driven models, and more.

Laboratory of Applied Mathematical Programming and Statistics (LAMPS PUC-Rio)

Who am I?

Hey there! I’m Matheus Nogueira, a Brazilian master’s student in Operations Research at LAMPS, an Optimization and Statistics Lab at PUC-Rio. I did my undergrad in Computer Engineering, also at PUC-Rio, in the Informatics’ Department.

Most of my work for the last years is about making sense of time series — mainly trying to predict the future with forecasting and simulation and I would love to share some interesting tools and concepts with you.

What’s LAMPS?

LAMPS is a research lab that mixes optimization, statistics, and machine learning — with people from both Electrical and Industrial Engineering at PUC-Rio. Besides academic research, we also work with companies and develop open-source packages, mostly in Julia.

If you’re into time series, you’ll probably like our GitHub. We’ve got packages for ARIMA, state space models, score driven models, and more.

Laboratory of Applied Mathematical Programming and Statistics (LAMPS PUC-Rio)

Why this article?

State space models are awesome. They’re flexible, powerful, intuitive, and cover a huge range of time series models — from the good old ARIMA to structural models and dynamic regression with exogenous variables.

In this article, I’ll walk you through how to use StateSpaceModels.jl, a package we’ve built at LAMPS to make modeling and forecasting with state space models easy and fast.

Home · StateSpaceModels.jl

We’ll do a hands-on tutorial, with some theory here and there to help you understand what’s going on. By the end, you should be able to run your own forecasting project in Julia using `StateSpaceModels.jl` — from start to finish.

Feel free to follow along with the available code at my GitHub.

https://github.com/MathNog/MathNogMediumArticles

Let’s dive in 🚀

Setting things up

Before we jump into modeling, let’s get our Julia environment ready.

Here are the packages you’ll need:

• StateSpaceModels.jl — for building and estimating the models,
• DataFrames — to handle our time series data,
• CSV — to read datasets from files (even though we’ll load most stuff from built-in datasets),
• Plots — for, well… plotting stuff,
• Statistics.jl — for simple descriptive statictics (aka, taking the mean of a vector).

You can install everything by running:

import Pkg
Pkg.add("StateSpaceModels")
Pkg.add("DataFrames")
Pkg.add("CSV")
Pkg.add("Plots")
Pkg.add("Statistcs)

If you’re working inside a project (which is a good idea for reproducibility), don’t forget to activate and instantiate the environment:

Pkg.activate(".")
Pkg.instantiate()

Then, just load the packages…

using StateSpaceModels, DataFrames, CSV, Plots, Statistics

If you’re following along using the notebook version of this tutorial, we also have some helper functions in a separate file:

include("statespacemodels_functions.jl")

Loading AirLinePassengers Dataset

Before jumping into model fitting, let’s take a look at a classic time series: the number of monthly airline passengers in the UK from 1949 to 1960. Fortunately, the package already has some built-in datasets.

Let’s load the dataset:

airline = CSV.File(StateSpaceModels.AIR_PASSENGERS) |> DataFrame;
airline.passengers = Float64.(airline.passengers);

Visual inspection

This dataset is probably one of the first things you see in any time series course — and for a good reason. It’s clean, has a strong upward trend, obvious seasonality, and some subtle challenges that make it perfect for learning.

But we’re not just going to throw a SARIMA at it and call it a day. What makes StateSpaceModels.jl really shine is its ability to handle structural models, where we decompose the series into unobserved components like trend and seasonality.

This is perfect for building models that are interpretable and flexible — and that’s exactly what we’ll explore soon.

First, let’s take a look into the series.

title = "Airline Passengers"
x_label = "Month"
y_label = "Passengers"
test_length_airline = 12*3
plot_time_series(log_airline, "month", "passengers", title, x_label, y_label; test_length = test_length_airline)

When you plot the original series, you’ll notice that the seasonal ups and downs get bigger over time. This means the variance of the seasonal component isn’t constant — it grows over time. That’s a problem for many models that assume the variance is stable (a.k.a. homoscedastic).

To deal with this, we apply a simple but powerful trick: log-transform the data. Taking the logarithm helps stabilize the unconditional variance, which makes it easier for our models to learn the structure of the series — especially the seasonal part. If, however, we have evidence of conditional heteroscedasticity, the log itself is not sufficient to handle it.

Here’s how we load and log-transform the dataset and plot the transformed series. See how the ups and downs seems to have the same amplitude over time?

log_airline = deepcopy(airline);
log_airline.passengers = log.(log_airline.passengers);

title = "Airline Log-Passengers"
x_label = "Month"
y_label = "Log(Passengers)"
test_length_airline = 12*3
plot_time_series(airline, "month", "passengers", title, x_label, y_label; test_length = test_length_airline)

Train-test split

As in any good machine learning workflow, we shouldn’t train our model on all the data and then judge its performance on the same data. That’s cheating.

To really know how well our model is doing, we need to test it on data it hasn’t seen before — data that wasn’t used during estimation. That’s why we split our dataset into two parts:

- a training set, to fit the model,

- and a test set, to evaluate its forecasting accuracy.

In time series, this split has to respect the order of the data — no shuffling allowed! We simply cut the last N observations to serve as our test set and keep the rest for training. The images above have already shown the test set for our tutorial.

airline_train, airline_test = time_series_split(log_airline, test_length_airline);

A quick dive into structural time series models

At the heart of StateSpaceModels.jl is the state space representation of time series models — a flexible way to write models that evolve over time. Some classic time series models, such as the ARIMA family and Exponential Smoothing (ETS) family can be written in this form. The state space framework and the Kalman Filter as strongly related, but this is a topic for another time. For this tutorial, we can simply see the Kalman Filter as a black-box method that allows us to fit a model that has a state space form.

The only theory that I wish to present is the general state space representation as a system of observation and state equations. For any given model, if we can write it in this format, we can say that this model has a state space representation.

One particularly useful class is the structural time series models, where we explicitly model components like trend, seasonality, cycles, and noise.

The general form of a structural model is:

This is called the observation equation— it links the hidden structure of the model to the actual data we observe. In this case, the time series data is written as a sum of trend/level, seasonality, cycle and noise components.

Each of the components evolves over time according to its own state equation. For example, the Basic Structural Model (BSM), which we’ll use a lot here, is made of one observation equation and 3 state equations, for level, slope and seasonality. Given the BSM, the state would be a vector composed of the level, slope and seasonality components.

Basic Structural Model equations

The BSM state equations are composed of a random walk with stochastic drift (a.k.a local linear trend model) and a stochastic seasonality component.

This setup allows each component to evolve in time with its own randomness — making the model highly flexible.

When we fit a model, we obtain fitted values for each component variance, the observation equation variance, and each component itself.

The package supports several types of models that can all be written in state space form, including: SARIMA, ETS, DAR, BasicStructual, BasicStructuralExplanatory, LinearRegression and some naive models as well.

We’ll mostly focus on the BasicStructural model here — but once you get the hang of it, switching to other types is straightforward.

We have two ways to define a BSM model in StateSpaceModels.jl. First, with the UnobservedComponents class, where we can specify each component. Notice that, since our time series has a monthly frequency, it’s seasonality period is 12 (months).

ssm_airline = StateSpaceModels.UnobservedComponents(airline_train.passengers;
                                                    trend = "local linear trend",
                                                    seasonal = "stochastic 12",
                                                    cycle = "no");

We could also use the BasicStructural class.

seasonal_period = 12
ssm_airline = StateSpaceModels.BasicStructural(airline_train.passengers, seasonal_);

Fitting the model

Once we’ve defined our model, fitting it is super straightforward. The `fit!` function takes care of estimating the parameters — namely, the variances of the components: irregular, trend, slope, and seasonal as well as the initial states— using maximum likelihood and the Kalman.

If any of these variances are estimated as zero, it means the corresponding component is deterministic, not stochastic. In other words, the model believes that component doesn’t need randomness to explain the data.

Here’s how we do it:

StateSpaceModels.fit!(ssm_airline);

We can print the fitted results nicely with:

print_results(ssm_airline);

This shows us the estimated variances, standard errors, and p-values for each component. In the output above, we can see that both the slope and seasonal variances were estimated as zero. At first glance, this suggests these components might be deterministic. But let’s not rush — we should always look at the p-values of the related z-test that tests if each parameter is statistically equal to zero.

The p-values associated with the irregular, trend and seasonality variances are all smaller than 0.01, which means that the null-hypothesis of each z-test (which says that the given parameter is equal to zero) can be strongly rejected. That being said, even though the fitted variance of the seasonal component was zero, the t-test says that we have strong evidence that it is not, so that this component is actually stochastic. The slope component, on the other hand, seems deterministic, since the p-value from the t-test is 1, which means that we do not have evidence to conclude that it’s value is not zero.

Model diagnostics

A crucial step after fitting any model is to check the residuals — and in the state space framework, we do this by analyzing the standardized innovations.

These innovations are the fitted values of the error terms — especially the one from the observation equation. In other words, they represent what the model couldn’t explain, scaled by the uncertainty at each point.

To get them, we run the Kalman filter and extract the innovations:

filter_airline = StateSpaceModels.kalman_filter(ssm_airline)
std_inov = StateSpaceModels.get_standard_innovations(filter_airline);
std_residuals = StateSpaceModels.get_standard_residuals(filter_airline);

We can plot the standard innovations and the standard residuals to see that they are the same thing.

plot(airline_train.month, std_inov, title = "Standardized Innovations", label = "", xlabel = "Month", ylabel = "Standardized Innovations")

plot(airline_train.month, std_residuals, title = "Standardized Residuals", label = "", xlabel = "Month", ylabel = "Standardized Residuals")

Diagnostics Plot

The BSM model (and most state space models) assumes that the residuals are: uncorrelated over time and normally distributed (i.e., errors are Gaussian white noise).

We can easily check these assumptions with a built-in function called plotdiagnostics:

p_diag = StateSpaceModels.plotdiagnostics(kalman_filter(ssm_airline), size = (1000, 600), 
                                        plot_title = "BSM Airline Std Innovations Diagnostics",
                                        title = ["Residuals" "ACF" "Histogram" "Q-Q Plot"])

This will give us:

• A time series plot of the residuals,
• The autocorrelation function (ACF),
• A histogram to check symmetry and outliers,
• A Q-Q plot to visually assess normality

Getting the predictive, filtered and smoothed states

One of the best things about using structural state space models is that we don’t just get a black-box prediction — we actually see the inner structure of the time series.

In particular, we can access the latent states estimated by the Kalman filter, which give us the values (and uncertainty) of each hidden component: trend, slope, seasonality, etc.

But first: what exactly is a state?

In this context, the state is the collection of unobserved variables that drive the dynamics of the series. For the Basic Structural Model, that includes the level, slope, and seasonal effects.

We can estimate these states in three different ways:

- Predictive states: Estimated before seeing the current data point:

- Filtered states: Estimated after seeing the current data point

- Smoothed states: Estimated using all the data (past and future)

You can get each one with:

predictive_states = StateSpaceModels.get_predictive_state(ssm_airline);
filtered_states   = StateSpaceModels.get_filtered_state(ssm_airline);
smoothed_states   = StateSpaceModels.get_smoothed_state(ssm_airline);

From the estimated states, we can extract the actual time series components — like the smoothed trend, filtered slope, or predictive seasonality:

smoothed_level, smoothed_slope, smoothed_seasonality       = get_components_from_states(smoothed_states);
filtered_level, filtered_slope, filtered_seasonality       = get_components_from_states(filtered_states);
predictive_level, predictive_slope, predictive_seasonality = get_components_from_states(predictive_states);

Why is this important?

Because now we can study the time series itself, not just predict it with a black-box model (yes, neural networks, I am talking to you).

Using the filtered and smoothed components, we can answer questions like:

• How did the trend evolve over time?
• Was there a sudden shift in growth?
• Did seasonality change or stay stable?
• When did the series deviate from its usual behavior?

These insights are especially useful in fields like economics, energy, retail, or any context where understanding the structure is as important as forecasting the future.

There is a built-in method to plot the smoothed components, which can be used as follows:

ks = StateSpaceModels.kalman_smoother(ssm_airline)
plot(ssm_airline, ks)

If we want to see the other components, we should do it by hand, as it is done below.

A curious mind may be asking itself why the initial values from the slope component seems unstable in the beginning of the series. This usually happens with models fitted with the Kalman filter, which needs the first observations to “warm up”. This is so common that in the following functions, I have chosen to exclude the first 12 observations from all components. This is not cheating if you know what you are doing!

plot_components(airline_train.month, predictive_level, predictive_slope, predictive_seasonality, "Predictive")

p_filtered = plot_components(airline_train.month, filtered_level, filtered_slope, filtered_seasonality, "Filtered")

It would be important to guarantee that the output we saw earlier, with the fitted variances, correspond to nature of the plotted components. Since the fitted BSM is a predictive model, we must analyse the predictive components. We have some intresting results! First, it is clear that the level/trend component is stochastic. The slope component is tricky. We expected a deterministic component, since it's variance was estimated as zero. However, the plot itself does no show this result and one could argue that it is stochastic. However, I would argue that even though the component is changing over time, is scale, in comparison to the series itself, is irrelevant! Tha is an interesting discussion to show that interpreting a model is not black and white, there are grey areas! The seasonal component, although was fitted as stochastic, is almost deterministic!

Out-of-sample forecast and simulation

Alright, now for the main event 🥁
As a forecaster guy™, this is where things get exciting: how do we make forecasts from a fitted state space model?

Not just point forecasts — we also want probabilistic scenarios, because in real life, we care about uncertainty too.

Luckily, StateSpaceModels.jl makes this super easy. You only need three functions:

- forecast: gives you the distribution of the future states, namely the expected value and the covariances.
- simulate_scenarios: samples possible future paths from the fitted predictive distribution.
- forecast_expected_value: gives you the expected value (mean) of the forecast.

Let’s run them:

S = 1000
scenarios_airline      = simulate_scenarios(ssm_airline, test_length_airline, S);
forecast_airline       = forecast(ssm_airline, test_length_airline);
expected_value_airline = forecast_expected_value(forecast_airline);

Now you have:

• forecast_airline: full probabilistic forecast object (including variances).
• scenarios_airline: a bunch of random paths your series might follow.
• expected_value_airline: the most likely forecast path (mean prediction).

This is where structural models shine: they don’t just spit out a number — they give you a full distribution, respecting all the structure you’ve captured (trend, seasonality, uncertainty, etc.).

Next up, let’s visualize those forecasts and compare them to the actual test data with the following function (available in my github). It seems quite accurate doesn't it?

plot_forecast(airline_train, airline_test, vec(expected_value_airline), scenarios_airline[:, 1, :], "Airline Passengers", "month", "passengers", "Month", "Passengers");

An astute observer may notice that this is not the original series. What use do we have for the expected log-passengers values? None! We should always recover the original series scale after forecasting. However, we wont do it here…

Model comparison

Since there is a number of available models in the package, let's fit and forecast an AutoArima and AutoETS models (yes, there is even automatic forecasting modeling in the package!) to compare their forecasting accuraccy.

They are quite simple to use.

First, we need to define and fit the models (in this case there is no need to explicitly call the fit! method).

auto_arima_airline = StateSpaceModels.auto_arima(airline_train.passengers; seasonal = 12);
auto_ets_airline   = StateSpaceModels.auto_ets(airline_train.passengers; seasonal = 12);

After that, we can compute the forecasts and simulations.

scenarios_arima_airline       = simulate_scenarios(auto_arima_airline, test_length_airline, S);
expected_value_arima_airline  = forecast_expected_value(forecast(auto_arima_airline, test_length_airline));

scenarios_ets_airline        = simulate_scenarios(auto_ets_airline, test_length_airline, S);
expected_value_ets_airline   = forecast_expected_value(forecast(auto_ets_airline, test_length_airline));

And compute the mean MAE and CRPS across the full horizon.

mean_bsm_crps = mean(StateSpaceModels.evaluate_crps(airline_test.passengers, scenarios_airline[:, 1, :]));
mean_bsm_mae  = mean(StateSpaceModels.evaluate_abs_error(airline_test.passengers, vec(expected_value_airline)));

mean_arima_crps = mean(StateSpaceModels.evaluate_crps(airline_test.passengers, scenarios_arima_airline[:, 1, :]));
mean_arima_mae  = mean(StateSpaceModels.evaluate_abs_error(airline_test.passengers, vec(expected_value_arima_airline)));

mean_ets_crps = mean(StateSpaceModels.evaluate_crps(airline_test.passengers, scenarios_ets_airline[:, 1, :]));
mean_ets_mae  = mean(StateSpaceModels.evaluate_abs_error(airline_test.passengers, vec(expected_value_ets_airline)));

DataFrame(
    Model = ["BSM", "ARIMA", "ETS"],
    CRPS = [mean_bsm_crps, mean_arima_crps, mean_ets_crps],
    MAE = [mean_bsm_mae, mean_arima_mae, mean_ets_mae]
)

We can see that, for this test-set, the BSM was the worst one with respect to the probabilistic metric (CRPS), while being the second best with respect to the point forecast metric (MAE).

Cross-validation

Let’s wrap things up with a super useful tool: cross-validation.

If you come from machine learning, you probably know this idea already — we split the data into multiple training and testing segments, fit the model on each training chunk, and evaluate how well it performs on the corresponding test chunk.

But in time series, we can’t just randomly shuffle data around. Time matters.

That’s why we use rolling or expanding window cross-validation. We train the model on a moving window of past data, then test it on a fixed-size future window. Repeat this several times and average the results.

This gives you a much more robust estimate of forecasting performance — especially when your series changes over time.

The image below, taken from another Medium article, illustrates a similar concept, in which each window incorporates new observations, namely the expanding window cross validation.

Cross Validation in Time Series

In StateSpaceModels.jl, it’s as easy as:

start_idx = 50
cv = StateSpaceModels.cross_validation(ssm_airline, test_length_airline, start_idx);

where start_idx represents the first observation of the first window to fit the model.

This function returns mean MAE and CRPS metrics from all windows by step ahead. This is an interesting plot because it enables us to see how the model forecast performance evolve along the forecast horizon. As we can see below, the forecast accuracy decreases over time. Fair enought! Further in the future, more difficult to foreacst!

plot_cv = plot(cv, "Basic Stuctural Model", size = (1000, 600))

Thanks for reading this far! I hope this tutorial helped you understand how powerful and intuitive `StateSpaceModels.jl` can be when it comes to modeling, interpreting, and forecasting time series in Julia. We went from loading data to fitting structural models, diagnosing residuals, extracting components, forecasting future values, and even cross-validating — all within one clean and consistent framework. There are even mode features available in the package that must be checked out!

Please refer to the article in Arxiv for more information regarding the package.

StateSpaceModels.jl: a Julia Package for Time-Series Analysis in a State-Space Framework

All the code used here is 100% written by me and available on my GitHub. The text, on the other hand, was written with the help of AI — because good tooling helps with writing too 😉

If you liked this package, be sure to check out the other time series tools we’ve built at LAMPSPUC. And stay tuned — more tutorials are on the way!

https://github.com/MathNog/MathNogMediumArticles

Laboratory of Applied Mathematical Programming and Statistics (LAMPS PUC-Rio)