Fitting a line: least squares

You have a cloud of noisy data points and you suspect a straight-line trend. Linear regression finds the one line y = m x + b that fits best — and "best" has a precise meaning: the line that makes the sum of squared vertical distances from the points as small as possible.

Remarkably, that best line has a simple closed-form formula (no searching required). Below you can set the underlying trend and how much noise we sprinkle on, then watch the fitted line track the data live in WebAssembly, alongside the true line it's trying to recover.

begin
    using PlutoUI, WasmMakie
end

true slope = 15 ÷10

noise level = 12 ÷10

number of points = 45

seed = 5

let
    # build noisy data y = m x + b + noise, accumulating the sums needed for the
    # closed-form least-squares fit in the SAME single pass
    m_true = Float64(slopei) / 10.0
    b_true = 2.0
    noise = Float64(noisei) / 10.0
    n = npoints
    s = (seed * 2654435761 + 12345) % 2147483647
    if s == 0
        s = 1
    end
    dx = Vector{Float64}(undef, n)
    dy = Vector{Float64}(undef, n)
    sx = 0.0
    sy = 0.0
    sxx = 0.0
    sxy = 0.0
    for i in 1:n
        x = 10.0 * Float64(i - 1) / Float64(n - 1)
        s = (s * 16807) % 2147483647
        u = Float64(s) / 2147483647.0
        y = m_true * x + b_true + noise * (u - 0.5) * 2.0
        dx[i] = x
        dy[i] = y
        sx = sx + x
        sy = sy + y
        sxx = sxx + x * x
        sxy = sxy + x * y
    end
    nn = Float64(n)
    m_fit = (nn * sxy - sx * sy) / (nn * sxx - sx * sx)
    b_fit = (sy - m_fit * sx) / nn

    fig = Figure(size = (600, 360))
    ax = Axis(fig[1, 1])
    # the data points, each drawn as a small vertical tick
    for i in 1:n
        lines!(ax, [dx[i], dx[i]], [dy[i] - 0.25, dy[i] + 0.25])
    end
    # the true line (what generated the data) and the fitted line (what we recovered)
    lines!(ax, [0.0, 10.0], [b_true, m_true * 10.0 + b_true])
    lines!(ax, [0.0, 10.0], [b_fit, m_fit * 10.0 + b_fit])
    fig
end

fit_stats = let
    # compute the least-squares fit in a SEPARATE bond-dependent cell so the markdown
    # below interpolates it live (values inside a markdown cell's own `let` bake to the
    # slider defaults).
    m_true = Float64(slopei) / 10.0
    noise = Float64(noisei) / 10.0
    n = npoints
    s = (seed * 2654435761 + 12345) % 2147483647
    if s == 0
        s = 1
    end
    sx = 0.0
    sy = 0.0
    sxx = 0.0
    sxy = 0.0
    for i in 1:n
        x = 10.0 * Float64(i - 1) / Float64(n - 1)
        s = (s * 16807) % 2147483647
        u = Float64(s) / 2147483647.0
        y = m_true * x + 2.0 + noise * (u - 0.5) * 2.0
        sx = sx + x
        sy = sy + y
        sxx = sxx + x * x
        sxy = sxy + x * y
    end
    nn = Float64(n)
    m_fit = (nn * sxy - sx * sy) / (nn * sxx - sx * sx)
    b_fit = (sy - m_fit * sx) / nn
    (floor(m_fit * 1000.0) / 1000.0, floor(b_fit * 1000.0) / 1000.0, floor(m_true * 100.0) / 100.0)
end;

Fitted line: y = 1.486 x + 2.125 (true slope 1.5, true intercept 2.0). Add more noise and the fit wobbles; add more points and it locks back onto the truth – more data beats noisier data.

Where the formula comes from

"Best fit" means minimizing the total squared error Σ (yᵢ − m xᵢ − b)². Setting the derivatives with respect to m and b to zero gives two linear equations — the normal equations — whose solution is exactly the slope and intercept we computed above. No iteration, no learning rate: a single formula.

This is the simplest member of a huge family. Replace the line with a plane, a polynomial, or millions of features and the same least-squares idea becomes the backbone of statistics and machine learning. It is also the moment "data science" stops being plotting and starts being modeling: proposing a relationship and measuring how well it explains the data.

Appendix

The MIT lecture fits with GLM.jl over DataFrames. WebAssembly can't run that stack, so we generate the data with an inline Park-Miller generator and apply the closed-form least-squares formula by hand, drawing with WasmMakie. The fit is identical to what GLM would return.