The chapter that closes the course is the one that holds the whole course together. Every model seen so far — the ARIMA coefficients of chapter 4, the weights of the MLP of chapter 5, the hyperparameters of the tree ensembles of chapter 6 — is a set of parameters, and the act of choosing them is parameter fitting. The deck opens with the problem statement in one sentence: we have a series of data points and a parametric mathematical function; we want to find the parameter values that maximize the fit of the function to the points.
The topic, the deck warns, is old and crowded: it has been studied, largely independently, in several disciplines — mathematics, statistics, computer science — and it goes under different names: parameter estimation, parameter fitting, curve fitting, (regression analysis). The formal umbrella is estimation theory, "a branch of statistics that deals with estimating the values of parameters based on measured empirical data that has a random component… An estimator attempts to approximate the unknown parameters using the measurements". The deck's promise: we will fly low, only scratching some theory.
| Name | Home discipline | Same act, different emphasis |
|---|---|---|
| Parameter estimation | statistics | the parameters describe an underlying physical setting; the estimator is judged on bias, variance, consistency |
| Curve fitting | numerical analysis | the geometric act: make a curve pass close to the points |
| Regression analysis | statistics / econometrics | the relationship between a response and predictors, with the noise modelled |
| Training / learning | machine learning | the same minimization, on huge parameter spaces, with gradients |
And then the deck says the sentence that positions this chapter in Part IV: as it is often the case, at the core there is optimization. We want to minimize the difference between the sample points and those predicted by the model function; the model function is given, we want to minimize its prediction errors. Formally: we want to minimize an error (loss, cost) function whose domain is the parameter space and whose image quantifies the quality of the prediction. Everything that follows — analytical least squares, gradient descent, simplex search, swarms, Optuna — is a different way of doing that one minimization.
Every forecasting chapter already did parameter fitting in disguise: the loss functions of chapters 2–7 are error functions over parameter spaces, and backpropagation is gradient descent on a very large parameter space. This chapter says it out loud and adds the tools for the cases where no gradient exists. The exam project is, at its core, one big fitting exercise — which is why the deck's examples (the Gaussian, the Box-Jenkins series) are the same ones used in the exam.
Given sample data (x₁, …, xₙ) and model data (ŷ ₁, …, ŷₙ), the quality of a fit is measured by an error function. The deck lists the standard family — BIAS, the arithmetic average of the errors; MAD, the Mean Absolute Deviation; MSE, the Mean Squared Error; the Standard Error; and MAPE, the Mean Absolute Percent Error.
The formulas on the deck's error-functions slide are printed images that did not survive text extraction. The definitions below follow the deck's own forecast_accuracy code, which is the authoritative statement of what each name means in this course.
| Metric | Definition | Notes |
|---|---|---|
| BIAS (ME) | ∑(ŷ − y) / n | mean error; can be ~0 on visibly wrong fits |
| MAD (MAE) | ∑|ŷ − y| / n | mean absolute deviation; robust to outliers |
| MSE | ∑(ŷ − y)² / n | mean squared error; punishes large errors quadratically |
| RMSE | √MSE | same unit as the data; the workhorse of the course |
| Standard error | standard deviation of the errors | √(MSE − BIAS²); error spread around its own mean |
| MAPE | mean(|ŷ − y| / |y|) × 100 | scale-free percentage; explodes when y ≈ 0 |
import os, numpy as np, pandas as pd, matplotlib.pyplot as plt
# from statsmodels.tsa.stattools import acf
# Accuracy metrics
def forecast_accuracy(forecast, actual):
mape = np.mean(np.abs(forecast - actual)/np.abs(actual)) # MAPE
me = np.mean(forecast - actual) # ME
mae = np.mean(np.abs(forecast - actual)) # MAE
mpe = np.mean((forecast - actual)/actual) # MPE
rmse = np.mean((forecast - actual)**2)**.5 # RMSE
corr = np.corrcoef(forecast, actual)[0,1] # corr
mins = np.amin(np.hstack([forecast[:,None], actual[:,None]]), axis=1)
maxs = np.amax(np.hstack([forecast[:,None], actual[:,None]]), axis=1)
minmax = 1 - np.mean(mins/maxs) # minmax
# acf1 = acf(forecast-actual)[1] # ACF1
return({'mape':mape, 'me':me, 'mae': mae, 'mpe': mpe,
'rmse':rmse, 'corr':corr, 'minmax':minmax})
if __name__ == '__main__':
os.chdir(os.path.dirname(os.path.abspath(__file__)))
df = pd.read_csv('whichever.csv') # can try with testLinear.csv
print( forecast_accuracy(forecast_values, actual.values) )Eleven points generated from y = 3 + 0.6x plus noise. Move the two coefficients and watch each metric react: the residual lines are the deviations the error functions average.
Which error function you report changes what "best" means. MSE/RMSE are the ones the least squares machinery minimizes — the square makes one large error dominate ten small ones, so the optimum chases the outliers. MAD is robust to them. MAPE is scale-free, which makes it portable across series of different magnitudes (chapters 2–7 used it for exactly that reason), but it is undefined when an actual value is zero. State your metric when you compare models in the exam project; the comparison of chapter 8 needs a defined criterion, and the Diebold–Mariano test takes it as an argument (crit="MSE").
The simplest case is also the founding one. We have a set of data points in ℝn (on this slide in ℝ2), we assume a linear model ŷ = θ₀ + θ₁x₁ + … + θₙxₙ, and we choose as error function MSE — actually the sum of squared errors, the "no mean" version. The error, plotted against the coefficients, defines the error surface to minimize: a bowl whose bottom is the best fit.
The optimal set of coefficients is determined by minimizing the loss (cost, objective) function, which is the maximum likelihood estimator for the task of linear regression. For the 1D case — a single array of n points, looking for the interpolation line ŷ = θ₀ + θ₁x that minimizes the mean squared error over all data pairs (xᵢ, yᵢ) — the solution can be written in closed form. With ̅ the average of the x's and ŷ̄ the average of the y's:
θ₁ = ∑(xᵢ − x̄)(yᵢ − ȳ) / ∑(xᵢ − x̄)² (equivalently: (n∑xy − ∑x∑y) / (n∑x² − (∑x)²))
θ₀ = ȳ − θ₁ · x̄th1 = ((X*y).mean() - X.mean()*y.mean())/((X**2).mean()-(X.mean())**2)
th0 = y.mean() - th1*X.mean()Then the deck draws the line that motivates the rest of the chapter: if the number of elements in X increases, so does the load on CPU/GPU, maybe too much. We can approximate the solution by minimizing the mean squared error by gradient descent (GD). The closed form is exact but global — it touches every data point at once; gradient descent is iterative and local, and scales to problems where the closed form does not exist at all.
"Maximum likelihood" is the bridge to chapter 8: if the noise around the true line is Gaussian, the probability of observing the data is maximized exactly when the sum of squared errors is minimized. Least squares is not an arbitrary choice of error function — it is the fitting rule that the Gaussian noise model dictates. That single equivalence is why the Gaussian appears at both ends of this course: as the distribution of measurement error (ch. 8) and as the justification of the squared loss (here).
Before any optimizer, the deck spends three slides on the calculus it will lean on. Let f(x) be a real-valued, twice differentiable function. The value of x for which the first derivative f′(x) is 0 corresponds to a maximum, a minimum or possibly a flex of f(x): for a maximum the second derivative f″(x) is negative, for a minimum it is positive, and it is 0 at an inflexion point. For functions of many variables, the same ideas are carried by the Jacobian (all first partial derivatives) and the Hessian (all second partial derivatives). The tabs below are the deck's three slides in reference form.
The value of x for which the first derivative f′(x) is 0 corresponds to a maximum, a minimum or possibly a flex of f(x).
Throughout the deck, first-order derivatives of f(x) are denoted f′(x), J(x) or ∇(x), "according to my whims".
the running example: f(x) = x³ − 4x² + 6 f′(x) = 3x² − 8x
stationary points at x = 0 (maximum, f″ = −8) and x = 8/3 (minimum, f″ = +8)The Jacobian operator is a generalization of the derivative operator to vector-valued functions: the Jacobian matrix of a vector-valued function f is the matrix of all its first-order partial derivatives. The Jacobian matrix represents the differential of f at every point where f is differentiable.
J = [ ∂f₁/∂x₁ ∂f₁/∂x₂ ] f: ℝⁿ → ℝⁿ ⇒ J is n × n
[ ∂f₂/∂x₁ ∂f₂/∂x₂ ]For a scalar error function E(θ) the Jacobian degenerates into the gradient row — the vector of partial derivatives that gradient descent will follow in section 5.
The Hessian matrix is a square matrix of second-order partial derivatives of a function f; it describes the local curvature of a function of many variables. The trace of the Hessian matrix is known as the Laplacian operator. In the deck, second-order derivatives are denoted f″(x), H(x) or ∇²(x).
H = [ ∂²f/∂x₁² ∂²f/∂x₁∂x₂ ]
[ ∂²f/∂x₂∂x₁ ∂²f/∂x₂² ]For the deck's two-variable example f(x₁, x₂) = 2x₁² + x₂² + x₁x₂ − 6x₁ − 5x₂, the Hessian is constant: H = [[4, 1], [1, 2]], positive definite — the surface is a convex bowl and every descent path ends in the same minimum (section 10's widget fits exactly this function).
Gradient (or steepest) descent (Cauchy, 1847) is the core of many learning algorithms (machine learning, deep learning, …). GD is a first-order iterative algorithm for finding the local minimum of a differentiable function f(x) — it is gradient ascent for the maximum. To find a local minimum it starts from a point x⁰ and takes steps proportional to the negative of the gradient (or of an approximate gradient) of the function at the current point; the length of the steps is dictated by a parameter, the learning rate (λ):
xₘ₊₁ = xₘ − λ · f′(xₘ)The deck works the idea by hand on the cubic of the toolkit. Cost function f(x) = x³ − 4x² + 6, derivative f′(x) = 3x² − 8x (the slope of f at any given x). Start with any value of x, say 0.5, and λ = 0.05. Iterate a number of times, always re-calculating x = x − λ · f′(x):
x = 0.5 − (−3.25 × 0.05) = 0.6625
x = 0.6625 + (3.983 × 0.05) = 0.86165 and so on, for a number of iterations,
or until a predefined value of precision is reachedThen the two-variable version, which shows what "step size" means once the surface has a shape. The function f(x₁, x₂) = 2x₁² + x₂² + x₁x₂ − 6x₁ − 5x₂ has partial derivatives ∂f/∂x₁ = 4x₁ + x₂ − 6 and ∂f/∂x₂ = 2x₂ + x₁ − 5, and a generic descent loop:
def f(x1,x2):
return 2*x1*x1 + x2*x2 + x1*x2 -6*x1 - 5*x2
def dfdx1(x1,x2):
return 4*x1 + x2 - 6
def dfdx2(x1,x2):
return 2*x2 + x1 - 5
def gradfun(currpoint):
x1 = currpoint[0]
x2 = currpoint[1]
return np.array([dfdx1(x1,x2),dfdx2(x1,x2)])
def gradient_descent(gradient, start, learn_rate, n_iter):
currpoint = start
for _ in range(n_iter):
diff = -learn_rate * gradient(currpoint)
currpoint += diff
return currpointf(x) = x³ − 4x² + 6 from x₀ = 0.5. The deck's trace with λ = 0.05 is x₁ = 0.6625, x₂ = 0.86165, …; push the learning rate up and watch the same recursion oscillate or explode.
The deck's worked example only shows the happy path. The same recursion with λ = 0.25 or 0.30 oscillates across the minimum and eventually flies away — the step is proportional to the slope, and near a steep wall a large λ overshoots into a steeper region still. The learning rate is itself a parameter that must be fitted (section 9), which is the recursion this whole chapter keeps running into.
Applied to linear regression, the machinery of section 5 becomes the batch algorithm that trains every linear model in this course. The cost function is
J(θ₀, θ₁) = (1/n) ∑ᵢ (θ₀ + θ₁xᵢ − yᵢ)²with partial derivatives
∂J/∂θ₀ = (2/n) ∑ᵢ (θ₀ + θ₁xᵢ − yᵢ) · 1
∂J/∂θ₁ = (2/n) ∑ᵢ (θ₀ + θ₁xᵢ − yᵢ) · xᵢMoving in the opposite direction of the gradient reduces the error; the deck notes that this needs X as a couple (1, x) — the intercept is absorbed by a column of ones, so the whole update is one matrix product:
def gradient_descent(X,y,theta,learning_rate=0.01,iterations=100):
m = len(y)
cost_history = np.zeros(iterations)
theta_history = np.zeros((iterations,2))
for it in range(iterations):
prediction = np.dot(X,theta)
theta = theta - (1/m)*learning_rate*( X.T.dot((prediction - y)))
theta_history[it,:] =theta.T
cost_history[it] = cal_cost(theta,X,y)
return theta, cost_history, theta_historyTwelve points from y = 2.5 + 0.7x plus noise. The cobalt line is the closed-form solution of section 3; the vermilion line is gradient descent from (0, 0) after the chosen number of steps. The lower band tracks the cost.
The closed form is one exact pass over the data — for linear models it is unbeatable. Gradient descent wins when the data no longer fits in one pass (the "CPU/GPU load" the deck mentions), when the model is nonlinear and no closed form exists (section 7), or when the "data" are millions of examples and a batch is only an approximation (stochastic variants). The two lines on the same plot are the same answer reached by different machinery — a useful sanity check whenever you implement one of them.
The nonlinear case the deck works in full is fitting a Gaussian curve. The Gaussian function g(x; a, μ, σ) = a · exp(−(x − μ)² / (2σ²)) has three parameters: the amplitude a (the peak height), the location μ (where the peak sits) and the scale σ (the spread). The deck's remark on the slide: no σ in the amplitude coefficient — when the Gaussian is written this way, a is the peak height independent of the width, which is exactly why the derivative with respect to a is the plain exponential.
The derivative formulas on the deck's slides 18–20 are images that did not survive text extraction. The standard derivatives are given below; they are the ones that make the deck's own gradient_descent code (slide 21) work, and the widget of this section reproduces the deck's numbers with them.
Derivatives of the Gaussian with respect to its parameters:
∂g/∂a = exp(−(x − μ)² / 2σ²)
∂g/∂μ = a · exp(−(x − μ)² / 2σ²) · (x − μ) / σ²
∂g/∂σ = a · exp(−(x − μ)² / 2σ²) · (x − μ)² / σ³The error is computed as (g(xᵢ; a, μ, σ) − yᵢ)², then averaged over all data points, and gradient descent acts on the parameters of the Gaussian function in order to minimize the error function. The gradients of the error function to include in the algorithm are the chain rule applied to each derivative: dE/dθ = mean(2 (g − y) · ∂g/∂θ) for each of θ = a, μ, σ. The deck's loop:
def gradient_descent(x, y, start, lmbda, niter):
currsol = start
diff = np.zeros(3)
for it in range(niter):
ampl = currsol[0]
mu = currsol[1]
sigma = currsol[2]
diff[0] = -lmbda * np.mean( dcostda(x, y, ampl, mu, sigma) )
diff[1] = -lmbda * np.mean( dcostdloc(x, y, ampl, mu, sigma) )
diff[2] = -lmbda * np.mean( dcostdscale(x, y, ampl, mu, sigma) )
currsol += diff
if it % 1000 == 0:
print("[{0}] {1}".format(it,currsol))
return currsolThirty-three points from g(x; a = 4, μ = 2, σ = 0.8) plus noise. Gradient descent starts from (a, μ, σ) = (2, 0, 2); the vermilion curve is where the three parameters have been dragged after the chosen number of steps.
Gaussian regression is the bridge between three tools the course has used separately: the Gaussian of chapter 8 (a distribution with parameters μ and σ), the fitting libraries Fitter and distfit (which try candidate distributions and return their parameters), and SciPy's curve_fit of section 8. Same objective — minimize the squared error against a parametric curve — with different optimizers behind the scenes.
SciPy provides the curve_fit() function for curve fitting via nonlinear least squares. The function takes x and y data as arguments, besides the name of the mapping function to use; curve_fit() returns the optimal values for the mapping function — the coefficient values — and also a covariance matrix for the estimated parameters. Once fit, the optimal parameters and the mapping function produce the output for any arbitrary input:
from scipy.optimize import curve_fit
# fit curve
popt, _ = curve_fit(objective, x_values, y_values)
# ------------------------------------ scipy
best_vals, covar = curve_fit(gaussian, x, y, p0=init_vals)
print('CF best_vals: {}'.format(best_vals))
yfit = gaussian(x,best_vals[0],best_vals[1],best_vals[2])The covariance matrix is the statistical payoff: its diagonal gives the uncertainty of each fitted parameter, which is the bridge to chapter 8 — a fitted coefficient is a point estimate, and its standard error is the square root of the corresponding diagonal entry.
curve_fit is a local method: behind the scenes it is a trust-region refinement (Levenberg–Marquardt) that needs a starting point, and p0 is how you provide it. For the Gaussian of section 7, a starting guess with σ ≈ 0 or μ far outside the data range can send the fit to a degenerate bell or to a different local minimum. On rough surfaces, the honest workflow is several restarts from different p0 values — the same lesson section 9 draws for search in general.
The first three sections fitted continuous parameters of known models with gradients. But algorithms are very sensitive to their control parameters (aka hyperparameters) — the λ of section 5 is the canonical example — and optimizing parameter setting is an active area of research, with many methods, from simple to very complicated. The simplest one: grid search: define an interval for each parameter and test all combinations (integer parameters). The deck's example is the ARIMA model of chapter 4, with parameters p, d, q over the intervals p ∈ [0, 3], d ∈ [0, 2], q ∈ [0, 2] — 4 × 3 × 3 = 36 configurations, evaluated exhaustively:
for (p=0;p<=3;p++)
for(d=0;d<=2;d++)
for(q=0;q<=2;q++)
{ performance = ARIMA(p,d,q);
if (performance > best_performance) … etcGrid search is exact but blind: it costs kⁿ evaluations for n parameters with k values each, and it has no idea where the interesting region is. For continuous parameters, it is possible to randomly generate values inside the interval — but this could lead to uneven sampling in multidimensional spaces, with clusters of near-duplicate configurations and empty voids elsewhere. A possible solution is Hammersley sampling: Hammersley is a representative of methods called "low discrepancy methods", with nontrivial mathematics involved; other possibilities include Sobol and Halton sequences (available in SciPy). They generate a sequence of n-dimensional points, each of which could correspond to a parameter setting.
The deck's printed output for one run — the first rows show the pattern: each column cycles through its own base sequence, so configurations spread instead of clustering:
| trial | θ₁ | θ₂ | θ₃ | θ₄ | θ₅ |
|---|---|---|---|---|---|
| 0 | 0.00000000 | 5.00000000 | 0.00000000 | 1.00000000 | 0.00000000 |
| 1 | 0.33333333 | 7.85714286 | 1.87500000 | 1.47619048 | 3.00000000 |
| 2 | 0.66666667 | 6.42857143 | 3.75000000 | 1.95238095 | 6.00000000 |
| 3 | 1.00000000 | 9.28571429 | 0.62500000 | 2.42857143 | 9.00000000 |
| 4 | 0.00000000 | 5.71428571 | 2.50000000 | 2.90476190 | 12.00000000 |
| 5 | 0.33333333 | 8.57142857 | 4.37500000 | 1.09523810 | 15.00000000 |
| 6 | 0.66666667 | 7.14285714 | 1.25000000 | 1.57142857 | 18.00000000 |
| 7 | 1.00000000 | 10.00000000 | 3.12500000 | 2.04761905 | 0.42857143 |
Grid search exhausts every combination; random search escapes the lattice but wastes evaluations on clusters; low-discrepancy sequences (Hammersley, Sobol, Halton) get the same coverage per evaluation as random search but without the voids. This is the whole history of hyperparameter search in one plate — and the direct ancestor of Optuna's smarter sampling in section 12.
Gradient descent needs a differentiable error function; many error surfaces in practice are not. Derivative-free methods, several of them, exist — Nelder–Mead simplex is one, essentially a local search. The Simplex (Nelder–Mead) method is simple to program, does not require derivative data, and has no general theoretical proof that it works — but lots of empirical evidence.
The mechanics: keep track of n + 1 points in n dimensions — the vertices of a simplex (a triangle in 2D, a tetrahedron in 3D, etc.). At each iteration the simplex can move, expand, or contract; it is sometimes known as the amoeba method, because the simplex "oozes" along the function. The basic operation is reflection: probe the point mirrored across the centroid of the best n vertices from the worst vertex. Then the decision tree:
The method is fairly efficient at each iteration (typically 1–2 function evaluations per iteration), but can take lots of iterations; it is somewhat flaky — sometimes it needs a restart after the simplex collapses on itself. The benefits balance the account: simple to implement, doesn't need a derivative, doesn't care about function smoothness. The deck's pseudocode, from Wikipedia:
We are trying to minimize the function f(x), x ∈ ℝⁿ; current test points x₁, …, xₙ₊₁
1. Order according to the values at the vertices, f(x₁) ≤ f(x₂) ≤ … ≤ f(xₙ₊₁).
Check whether the method should stop.
2. Calculate xo, the centroid of all points except xₙ₊₁.
3. Reflection: xr = xo + α(xo − xₙ₊₁), α > 0.
If f(x₁) ≤ f(xr) < f(xₙ), replace the worst point with xr and go to step 1.
4. Expansion: if f(xr) < f(x₁), compute xe = xo + γ(xr − xo), γ > 1.
If f(xe) < f(xr), replace the worst with xe and go to step 1; else keep xr.
5. Contraction: here f(xr) ≥ f(xₙ). Compute xc = xo + ρ(xₙ₊₁ − xo), 0 < ρ ≤ 0.5.
If f(xc) < f(xₙ₊₁), replace the worst with xc and go to step 1.
6. Shrink: replace all points except the best (x₁) with xᵢ = x₁ + σ(xᵢ − x₁) and go to step 1.
α, γ, ρ and σ are the reflection, expansion, contraction and shrink coefficients; standard values α = 1, γ = 2, ρ = 1/2, σ = 1/2.The deck's own function f(x₁, x₂) = 2x₁² + x₂² + x₁x₂ − 6x₁ − 5x₂, minimum at (1, 2) with f = −8. The starting triangle is (5,5), (6,5), (5,6); the cobalt bands are the error surface, the triangle is the current simplex, and the label tells you which move the amoeba just performed.
Nelder–Mead is the local-search workhorse you will meet again in the prescriptive chapter, where decision spaces replace parameter spaces: the same ideas — move, expand, contract, shrink — reappear in metaheuristics. In SciPy it is one line: minimize(fun, x0, method='Nelder-Mead'). When the exam project asks why a model's hyperparameters were chosen this way, "I searched the space with a derivative-free local search / with Optuna" is the honest, complete answer — with the caveat that, being local, it needs restarts on rough surfaces.
The last search metaphor of the deck is biological. PSO is an optimization method inspired by the simulation of behaviors of animals that move in flocks; it uses social interaction as a metaphor to design a problem solving algorithm. Proposed in 1995 by James Kennedy (social psychologist) and Russell Eberhart (electrical engineer), PSO algorithms are especially useful for parameter optimization in continuous, multi-dimensional search spaces.
The machinery: PSO uses agents (particles) which are members of a swarm that moves in the search space looking for an optimal solution; each particle is considered as a point in an n-dimensional space that regulates its "flight" according to its flight experience and that of the other particles. The swarm is made of particles in a multidimensional space, each with a position and a speed; the particles fly in this hyperspace (the search space) and can do two things: memorize their best position (best solution found) and get to know the global best, and/or the best in their vicinity. The members of a swarm communicate good positions to each other, and accordingly change their speed, and subsequently positions. How much influence a point has is determined by its fitness — a measure assigned to a potential solution, capturing how good it is compared to all other solution points. The metaphors the deck draws: evolution, with the idea of "survival of the fittest"; social behavior, through a "follow the local leader" effect; and emergent pattern formation.
Three bests organize the flight. Each particle keeps track of its coordinates associated with the best fitness it has achieved so far — pbest. Another "best" value is obtained by any particle in the neighbors of the particle — lbest. The best value obtained by the whole population is the global best — gbest. PSO search consists of iteratively changing the velocity of (accelerating) each particle toward its pbest and gbest (global version of PSO) or lbest locations (local version); acceleration is weighted by a random term, with separate random numbers for acceleration toward pbest and gbest/lbest.
Formally, each particle modifies its position on the basis of its current position, its own current speed, the distance between the current position and pbest, and the distance between the current position and gbest (and possibly lbest):
vᵢₘ₊₁ = w · vᵢₘ + C₁ · rand() · (pbestᵢ − xᵢₘ) + C₂ · rand() · (gbest − xᵢₘ)
xᵢₘ₊₁ = xᵢₘ + vᵢₘ₊₁
where vᵢₘ is the velocity of particle i at iteration k, w the inertia weight,
rand() uniform random numbers between 0 and 1, xᵢₘ the current position.The weights are often dynamically updated, from diversification to intensification as the search runs:
cᵢ = cₘₘₙx − [ (cₘₘₙx − cᵢₘₓ) × iter / maxiter ]The algorithm in seven steps, as the deck lists it: 1) initialize the population — locations and velocities; 2) evaluate the fitness of each individual particle (pbest); 3) keep track of the swarm / neighborhood highest fitness (gbest / lbest); 4) modify velocities based on pbest and gbest / lbest; 5) update the particle positions; 6) terminate if the condition is met; 7) go to step 2. In pseudocode:
1. Randomly initialize particle positions x[i,j] and velocities v[i,j]
2. Initialize the global and local fitness to the worst possible
3. Loop until end_condition
4. evaluate the fitness f[i] of each particle
5. update the particle bests pbest[i,j] of each particle i
6. update the global (neighborhood) best gbest (lbest)
7. // update of particle velocity and position
for i = 1 to number of particles n do
for j = 1 to number of dimensions m do
R1 = uniform random number
R2 = uniform random number
v[i,j] = w*v[i,j] + C1*R1*(pbest[i,j] - x[i,j]) + C2*R2*(gbest[j] - x[i,j])
x[i,j] = x[i,j] + v[i,j]| Parameter | Meaning | Typical setting (deck) |
|---|---|---|
| number of particles | swarm size | usually between 10 and 50 |
| C₁ | importance of the personal best (pbest) | often C₁ + C₂ = 4 (empirical) |
| C₂ | importance of the global / neighborhood best (gbest / lbest) | same sum rule |
| velocity | step magnitude per iteration | too low → too slow; too high → too unstable |
Contrast the three optimizers seen so far on one axis: gradient descent exploits one point with exact local information; Nelder–Mead exploits a small population of n + 1 points with function values only; PSO explores with a whole population and random pulls. Exploration is the price you pay for robustness — PSO and its cousins are the methods that keep working when the error surface is rough, multimodal and derivative-free, which is precisely the situation the exam project's hyperparameter search is in.
The deck closes with the practical answer to the whole chapter: Optuna is a hyperparameter optimization framework designed for efficient, automated search over parameter spaces. It optimizes an objective function minθ f(θ) where θ are hyperparameters and f(θ) is the cost (loss) function to minimize; Optuna treats this as a black-box optimization problem. It defines a Study as the optimization process where all trials and results are stored, and Trials as single evaluations of the objective function, corresponding to single hyperparameter configurations:
study = optuna.create_study(direction="minimize")
def objective(trial):
x = trial.suggest_float("x", -10, 10)
return x**2The search space is defined implicitly inside the objective via the suggest_* methods — one call per parameter, which both samples a candidate value and records the parameter for the history:
| Method | Use | Example from the deck |
|---|---|---|
trial.suggest_float | continuous parameters, optionally on a log scale | suggest_float("learning_rate", 1e-3, 0.2, log=True) |
trial.suggest_int | integer parameters | suggest_int("max_depth", 2, 8) |
trial.suggest_categorical | discrete choices | suggest_categorical("booster", ["gbtree", "dart"]) |
Search then loops over parameters sampling (different strategies, default TPE), evaluating the objective and updating the internal model: study.optimize(objective, n_trials=100). Search can be pruned: trial.report(value, step) and trial.should_prune() — trials that are already hopeless are killed early, which is what makes Optuna cheap on expensive objectives like training an XGBoost model. The deck's full example fits exactly the workflow of the exam project: forecast the Box-Jenkins airline passengers series with XGBoost, tuning the hyperparameters with Optuna.
import optuna
import xgboost as xgb
import numpy as np, pandas as pd
from sklearn.metrics import mean_squared_error
import matplotlib.pyplot as plt
def create_dataset(series, n_lags):
X, target = [], []
for t in range(n_lags, len(series)):
X.append(series[t - n_lags:t])
target.append(series[t])
return np.array(X), np.array(target)
y = pd.read_csv("BoxJenkins.csv",usecols=["Passengers"])
log_y = np.log(y.values.flatten())
nlags = 12 # use past 12 months to predict next month
X_all, y_all = create_dataset(log_y, nlags)
split_idx = int(0.8 * len(X_all))
X_train, X_valid = X_all[:split_idx], X_all[split_idx:]
y_train, y_valid = y_all[:split_idx], y_all[split_idx:]
# Optuna objective
def objective(trial):
params = {
"objective": "reg:squarederror",
"eval_metric": "rmse",
"booster": "gbtree",
"n_estimators": trial.suggest_int("n_estimators", 100, 800),
"max_depth": trial.suggest_int("max_depth", 2, 8),
"learning_rate": trial.suggest_float("learning_rate", 1e-3, 0.2, log=True),
"min_child_weight": trial.suggest_int("min_child_weight", 1, 10),
"subsample": trial.suggest_float("subsample", 0.6, 1.0),
"colsample_bytree": trial.suggest_float("colsample_bytree", 0.6, 1.0),
"gamma": trial.suggest_float("gamma", 1e-8, 5.0, log=True),
"reg_alpha": trial.suggest_float("reg_alpha", 1e-8, 5.0, log=True),
"reg_lambda": trial.suggest_float("reg_lambda", 1e-8, 5.0, log=True),
"random_state": 666,
}
model = xgb.XGBRegressor(**params)
model.fit(X_train, y_train, eval_set=[(X_valid, y_valid)], verbose=False)
pred_valid = model.predict(X_valid)
rmse = np.sqrt(mean_squared_error(y_valid, pred_valid))
return rmse# Run Optuna
study = optuna.create_study(direction="minimize")
study.optimize(objective, n_trials=50)
print("Best RMSE:", study.best_value)
print("Best params:")
for k, v in study.best_params.items():
print(f" {k}: {v}")
# Fit the final model on train + valid
best_model = xgb.XGBRegressor(
objective="reg:squarederror",
eval_metric="rmse",
booster="gbtree",
random_state=42,
**study.best_params
)
best_model.fit(X_all, y_all, verbose=False)# Sliding forecast
def sliding_forecast(model, last_window, horizon):
window = last_window.copy()
preds = []
for _ in range(horizon):
next_pred = model.predict(window.reshape(1, -1))[0]
preds.append(next_pred)
window = np.roll(window, -1)
window[-1] = next_pred
return np.array(preds)
H = 12 # forecast next 12 months
last_window = log_y[-nlags:]
forecast_log = sliding_forecast(best_model, last_window, H)
# Undo the log transform
forecast = np.exp(forecast_log)
print(f"Forecast: {np.round(forecast, 1)}")
time = np.arange(len(y))
future_time = np.arange(len(y), len(y) + H)
plt.figure(figsize=(10, 5))
plt.plot(time, y, label="Observed")
plt.plot(future_time, forecast, label="Forecast")
plt.axvline(len(y) - 1, linestyle="--")
plt.legend()
plt.title("Airline Passengers: XGBoost + Optuna")
plt.show()The deck's printed listing has two slips that would crash the copy-paste: the final print statement uses a broken quote (f“Forecast…) and the sliding-forecast call refers to N_LAGS while the variable is named nlags. Both are corrected above; the flow is otherwise exactly the deck's.
This is the recommended workflow for the exam project's tuning step: define the objective as a validation error (RMSE on the held-out span), declare the space with suggest_* calls, run 50–100 trials, refit the best configuration on train + valid, and report study.best_value and study.best_params. It is the modern, defensible successor of every method in this chapter: it keeps grid search's bookkeeping, random search's flexibility, low-discrepancy sampling's coverage, and adds a model (TPE) that learns which regions of the parameter space are promising — plus pruning for expensive objectives.
We have a series of data points and a parametric mathematical function; we want to find the parameter values that maximize the fit of the function to the points. Since the model function is given, this reduces to minimizing an error (loss, cost) function whose domain is the parameter space and whose image quantifies the quality of the prediction: at the core there is optimization.
Parameter estimation (statistics), parameter fitting, curve fitting (numerical analysis) and regression analysis (statistics/econometrics). The formal umbrella is estimation theory: estimating the values of parameters based on measured empirical data that has a random component, where an estimator attempts to approximate the unknown parameters using the measurements. In machine learning the same act is called training.
BIAS (ME): mean of the errors ∑(ŷ − y)/n — can be ~0 on visibly wrong fits. MAD (MAE): mean of the absolute deviations — robust to outliers. MSE: mean of the squared errors — punishes large errors quadratically. RMSE: square root of MSE — same unit as the data. Standard error: standard deviation of the errors, √(MSE − BIAS²). MAPE: mean(|ŷ − y|/|y|) × 100 — scale-free, but explodes when an actual value is near zero.
If the noise around the true linear relationship is Gaussian, the probability of observing the data is maximized exactly when the sum of squared errors is minimized. Least squares is therefore not an arbitrary loss choice: it is the fitting rule that the Gaussian noise model dictates. The equivalence is why the Gaussian appears both as the measurement-error distribution (chapter 8) and as the justification of the squared loss.
θ₁ = ∑(xᵢ − x̄)(yᵢ − ȳ) / ∑(xᵢ − x̄)² (equivalently (n∑xy − ∑x∑y)/(n∑x² − (∑x)²)); θ₀ = ȳ − θ₁x̄. Python: th1 = ((X*y).mean() - X.mean()*y.mean())/((X**2).mean()-(X.mean())**2); th0 = y.mean() - th1*X.mean().
When the number of elements in X increases, the load on CPU/GPU grows, maybe too much. Gradient descent approximates the same minimum iteratively and scales to large data and to nonlinear models where no closed form exists. The two routes agree on the same line, which makes them a sanity check on each other.
The Jacobian is the generalization of the derivative to vector-valued functions: the matrix of all first-order partial derivatives, representing the differential of f. The Hessian is the square matrix of second-order partial derivatives, describing the local curvature of a function of many variables. The trace of the Hessian is the Laplacian operator. Deck notation: f′(x), J(x) or ∇(x) for first order; f″(x), H(x) or ∇²(x) for second order.
xₘ₊₁ = xₘ − λ · f′(xₘ): from a starting point x⁰, take steps proportional to the negative of the gradient, with the step length dictated by the learning rate λ. It is a first-order iterative algorithm for local minima (Cauchy, 1847); too small a λ is slow, too large a λ overshoots and diverges — the learning rate is itself a parameter that needs fitting.
f′(x) = 3x² − 8x. First step: f′(0.5) = −3.25, so x₁ = 0.5 − (−3.25 × 0.05) = 0.6625. Second: f′(0.6625) ≈ −3.983, so x₂ = 0.6625 + (3.983 × 0.05) = 0.86165. The minimum is at x = 8/3, where f ≈ −3.48; the recursion approaches it as long as λ is small enough.
With J(θ₀, θ₁) = (1/n)∑(θ₀ + θ₁xᵢ − yᵢ)²: ∂J/∂θ₀ = (2/n)∑(θ₀ + θ₁xᵢ − yᵢ) and ∂J/∂θ₁ = (2/n)∑(θ₀ + θ₁xᵢ − yᵢ)xᵢ. Moving in the opposite direction of the gradient reduces the error; in matrix form X carries a column of ones so the update is one product: θ − (1/m)λXᵠ(ŷ − y).
g(x; a, μ, σ) = a · exp(−(x − μ)²/(2σ²)) with amplitude a, location μ, scale σ — "no σ in the amplitude coefficient": a is the peak height. Derivatives: ∂g/∂a = exp(−(x − μ)²/2σ²); ∂g/∂μ = g · (x − μ)/σ²; ∂g/∂σ = g · (x − μ)²/σ³. The error gradients are dE/dθ = mean(2(g − y) · ∂g/∂θ) for θ = a, μ, σ.
curve_fit (SciPy) takes x and y data plus the name of the mapping function; it returns the optimal parameter values and a covariance matrix for the estimated parameters (whose diagonal gives each parameter's uncertainty). p0 is the starting point for the local nonlinear-least-squares refinement: on rough surfaces, different p0 can reach different local minima, so restarts are part of the honest workflow.
Grid search defines an interval for each parameter and tests all combinations — the deck's ARIMA example: p ∈ [0,3], d ∈ [0,2], q ∈ [0,2], 36 evaluations. It costs kⁿ evaluations for n parameters with k values each. For continuous parameters, random generation inside the intervals leads to uneven sampling in multidimensional spaces (clusters and voids). The remedy is low-discrepancy sampling: Hammersley, Sobol and Halton sequences (the latter two available in SciPy), which generate evenly spread n-dimensional points, each a candidate parameter setting.
It keeps track of n + 1 points in n dimensions — the vertices of a simplex (triangle in 2D). Each iteration: order the vertices by function value; compute the centroid of all except the worst; try reflection across the centroid; if the reflected point is the best so far try expansion; if reflection didn't help try contraction toward the centroid; if all else fails shrink the simplex around the best vertex. Coefficients: α = 1 (reflection), γ = 2 (expansion), ρ = 1/2 (contraction), σ = 1/2 (shrink). Typically 1–2 function evaluations per iteration; simple, derivative-free, smoothness-agnostic, but sometimes flaky (restart after collapse).
pbest: the best fitness a particle has achieved so far, with its coordinates; lbest: the best found in the particle's neighborhood; gbest: the best found by the whole population. Each particle accelerates toward pbest and gbest (or lbest) with random weights: vᵢₘ₊₁ = w·vᵢₘ + C₁·rand()·(pbestᵢ − xᵢₘ) + C₂·rand()·(gbest − xᵢₘ), then xᵢₘ₊₁ = xᵢₘ + vᵢₘ₊₁. Typical settings: 10–50 particles, C₁ + C₂ = 4 (empirical); velocity too low → too slow, too high → too unstable; weights often decay linearly from diversification to intensification.
A study is the whole optimization process where all trials and results are stored (optuna.create_study(direction="minimize")); a trial is a single evaluation of the objective, i.e. one hyperparameter configuration. The search space is declared implicitly inside the objective with suggest_float (optionally log=True), suggest_int and suggest_categorical. Sampling uses different strategies (default TPE, a model that learns which regions are promising). Pruning (trial.report(value, step) + trial.should_prune()) kills trials that are already hopeless, saving computation on expensive objectives.