Statistics & ML

Every statistics, machine-learning and estimation builtin that runs today. The list is read off the interpreter's own dispatch table, so nothing here is aspirational. Plain descriptive statistics — mean, std, var, median — are in Core Maths; this chapter starts at correlation and goes up.

Example data

A regression problem and a held-out set. Every example below runs against these.

X    = randn(120, 3, seed = 3)   # 120 samples, 3 features

beta = [2.0, -1.0, 0.5]
y    = X * beta + 0.1 * randn(120, seed = 4)
Xnew = randn(30, 3, seed = 5)
ynew = Xnew * beta

Descriptive statistics

Two conventions worth knowing before you compare a number against another tool: kurtosis follows MATLAB (normal = 3), not SciPy's default (fisher=True, normal = 0); and corr(x, y) is exactly cov(x, y) / (std(x) * std(y)), because all three share the N-1 denominator.

corr and cov both answer "how do these two vectors move together", just on different scales: cov keeps the original units (so it changes if you rescale a or b), while corr divides that out and always lands in [-1, 1], which is why corr is the one worth reporting and cov is the one worth using inside a formula. The example below builds b as a known mix of a and independent noise, prints both numbers, then plots the scatter so the number and the shape it describes sit side by side.

a = randn(200, seed = 11)
b = 0.8 * a + 0.6 * randn(200, seed = 12)
print("corr(a, b): {round(corr(a, b), 3)}")
print("cov(a, b):  {round(cov(a, b), 3)}")
scatter(a, b, color = "#4169e1", alpha = 0.6)
title("corr = {round(corr(a, b), 3)}, cov = {round(cov(a, b), 3)}")
xlabel("a")
ylabel("b")

corr(a, b): 0.691
cov(a, b):  0.515
· scatter: figure 1 (200 points)

mode, skewness and kurtosis describe the shape of one sample rather than a relationship between two. mode just answers "which value showed up most", useful on counts or categories where a mean would land between categories that do not average into anything meaningful. skewness and kurtosis are about the tails: a skewed sample has one tail longer than the other, and a sample with heavy or light tails relative to a normal distribution shows up as a kurtosis above or below 3. The example below takes nine scores with one obvious outlier and reports all three in one go:

scores = [3, 5, 5, 5, 6, 7, 7, 9, 20]
print("mode: {mode(scores)}")
print("skewness: {round(skewness(scores), 3)}")
print("kurtosis: {round(kurtosis(scores), 3)}")

mode: 5
skewness: 1.957
kurtosis: 5.693

The one outlier at 20 pulls the tail out to the right, which is exactly what a positive skewness means, and pushes kurtosis well past the normal distribution's 3.

Regression

ridge adds no intercept — append a column of ones yourself if you want one. ridge_model does add one, which is the difference between the two names.

polyfit and polyval are a matched pair: fit once, then evaluate the same coefficients anywhere you like, whether that is the original x values (to check the fit) or a finer grid (to draw a smooth curve). The example below fits a noisy cubic, prints the recovered coefficients, and overlays the fitted curve on the noisy points so you can see how well a degree-3 polynomial tracks a signal that really is degree 3:

xs = -3 to 3 step 0.25
truth = 0.5 * xs .^ 3 - 2 * xs
noisy = truth + 1.5 * randn(len(xs), seed = 7)
c = polyfit(xs, noisy, 3)
print("polyfit coefficients (highest degree first): {round(c, 3)}")
scatter(xs, noisy, color = "#94a3b8", label = "measured")
plot(xs, polyval(c, xs), color = "#4169e1", width = 2, label = "cubic fit")
legend()
title("polyfit degree 3")

polyfit coefficients (highest degree first): [0.475, 0.039, -1.982, -0.056]
· scatter: figure 2 (25 points)
· plot: figure 3 (25 points)

ridge is for when a plain OLS fit would be unstable — collinear columns, or more columns than rows — and a small amount of shrinkage toward zero buys back stability at the cost of a little bias. Because ridge never learns an intercept, the example below sets up a design matrix X with zero mean by construction (randn), so the "true" coefficients really are recoverable without one:

X    = randn(120, 3, seed = 3)
beta = [2.0, -1.0, 0.5]
y    = X * beta + 0.1 * randn(120, seed = 4)
c = ridge(X, y, 0.5)
print("ridge coef: {round(c, 3)}")   # close to [2.0, -1.0, 0.5]

ridge coef: [1.993, -1.002, 0.496]

Time series: autoregressive models

An AR(p) model forecasts its own future from its own fitted history, so predict takes a horizon rather than new data — the one place in this chapter where predict does not mean "apply to new rows".

Past the first p steps the recursion is conditioning on its own earlier forecasts, because an out-of-sample future is not available to condition on. That is not a shortcut — it is what a multi-step AR forecast is — but it does mean the uncertainty grows with the horizon in a way the returned numbers do not show.

x = zeros(500)
x[0] = randn(seed = 21)
x[1] = randn(seed = 22)
for t = 2 to 499
    x[t] = 0.6 * x[t - 1] - 0.3 * x[t - 2] + 0.1 * randn()
end for
m = ar_model(x, 2)
print("generating: [0.6, -0.3]")
print("recovered:  {round(m.coef, 3)}")

future = m.predict(30)
plot(460 to 499, x[460:499], color = "#0f172a", label = "measured")
plot(500 to 529, future, color = "#e11d48", width = 2, label = "forecast")
legend()
title("AR(2), 30 steps ahead")

generating: [0.6, -0.3]
recovered:  [0.489, -0.153]
· plot: figure 4 (40 points)
· plot: figure 5 (30 points)

How close the recovery gets depends on how much record you give it. At 500 samples it lands within a few hundredths; at 100 it does not, and the forecast inherits that error compounded once per step.

Clustering and dimensionality reduction

These describe the dataset they were given and have no .predict. Where a reusable model makes sense there is a _model sibling — see the fitted-model protocol.

Qu has no multiple return, so getting both labels and centroids means calling kmeans and kmeans_centers with the same seed. dbscan has no _model sibling on purpose: a point's cluster depends on the whole dataset's local density, so there is no honest .predict(Xnew) to write. eps is in X's own units and nothing is auto-scaled, which is the usual reason a first DBSCAN run returns one cluster or all noise.

G = vstack(randn(40, 2, seed = 31) + 4,
           randn(40, 2, seed = 32) - 4,
           randn(40, 2, seed = 33))
labels = kmeans(G, 3, seed = 5)
centers = kmeans_centers(G, 3, seed = 5)
xs = G[:, 0]
ys = G[:, 1]
for k in 0 to 2
    inside = labels == k
    print("cluster {k}: {len(xs[inside])} points")
    scatter(xs[inside], ys[inside], label = "cluster {k}")
end for
scatter(centers[:, 0], centers[:, 1], marker = "x", color = "#0f172a", size = 12, label = "centroids")
legend()
title("kmeans, k = 3")

cluster 0: 40 points
· scatter: figure 6 (40 points)
cluster 1: 39 points
· scatter: figure 7 (39 points)
cluster 2: 41 points
· scatter: figure 8 (41 points)
· scatter: figure 9 (3 points)

The three blobs really do have 40 points each; k-means recovers 40, 39 and 41, and the one it misplaces sits between two clusters. A mask indexes a vector directly -- xs[labels == k] -- which is how you draw one cluster at a time, since color = takes one colour for the whole series rather than one per point.

dbscan needs an eps sized to the data rather than a cluster count; pca/pca_components/pca_explained_variance describe the same fit three different ways; tsne is for a quick look at a handful of points, not a metric to report:

G = vstack(randn(40, 2, seed = 31) + 4,
           randn(40, 2, seed = 32) - 4,
           randn(40, 2, seed = 33))
d = dbscan(G, 1.5, 4)
print("dbscan labels: {distinct(d)}")

proj  = pca(G, 2)                    # the (n, 2) projected coordinates

comps = pca_components(G, 2)         # the (2, 2) loading matrix

ev    = pca_explained_variance(G, 2)
print("pca projected shape: {rows(proj)} x {cols(proj)}")
print("pca components:\n{round(comps, 3)}")
print("explained variance: {round(ev, 3)}")

small = vstack(randn(6, 2, seed = 41) + 3, randn(6, 2, seed = 42) - 3)
low = tsne(small, n_components = 2, perplexity = 3, seed = 1, max_iter = 100)
print("tsne output rows: {rows(low)}")

dbscan labels: (0, 1, 2)
pca projected shape: 120 x 2
pca components:
[0.694, 0.72; 0.72, -0.694]
explained variance: [0.969, 0.031]
tsne output rows: 12

G already carries three well-separated blobs, so eps = 1.5 finds them without also finding noise; a smaller eps on the same data returns mostly -1. pca on a two-column input is not reducing anything -- it is here to show the three functions agree with each other -- the real use is k < d.

The fitted-model protocol

Every estimator fits with a constructor call and returns a handle: named read-only fields, plus predict/score (and update/estimate for the filters) reached through ordinary method sugar.

m = ols_model(X, y)
print(round(m.coef, 3))            # a field, not a method call

yhat = m.predict(Xnew)
print(round(m.score(Xnew, ynew), 4))

[2, -1.005, 0.497]
1

predict, update, estimate, score, fit and pipeline are each one builtin that switches on the model's own .kind, not a separate builtin per model type. That is why the same predict means "apply to new rows" for a regression and "time update" for a Kalman filter.

Constructors

Four things about this table are easy to get wrong:

• kind= is never inferred. A 0/1 target and a small-integer count target look identical to the type system, so tree_model and friends make you say kind = "classification" when that is what you mean.
• gradient_boosting_model defaults to max_depth = 3, shallower than tree_model's 20. A boosting weak learner is supposed to be weak.
• svm_model is binary only. y must have exactly two distinct values; there is no one-vs-rest wrapper.
• One seed= reproduces a whole forest. It seeds a master stream that each tree draws its own seed from, so the forest is deterministic and no two trees share a stream.

Xc = randn(200, 2, seed = 41)
yc = zeros(200)
for i in 0 to 199
    if Xc[i, 0] + Xc[i, 1] > 0
        yc[i] = 1
    end if
end for
s = train_test_split(Xc, yc, test_size = 0.3, seed = 9)

function accuracy(model)
    yhat = model.predict(s.X_test)
    hit = yhat == s.y_test
    return len(yhat[hit]) / len(yhat)
end function

names = ("knn", "tree", "forest", "bayes")
scores = (
    accuracy(knn_model(s.X_train, s.y_train, 5, kind = "classification")),
    accuracy(tree_model(s.X_train, s.y_train, kind = "classification")),
    accuracy(random_forest_model(s.X_train, s.y_train, 40, seed = 3, kind = "classification")),
    accuracy(naive_bayes_model(s.X_train, s.y_train))
)
for i in 0 to 3
    print("{names[i]}: {round(scores[i], 3)}")
end for
bar(0 to 3, scores, color = "#4169e1", values = true)
xticks(0 to 3)
xticklabels(names)
ylim(0, 1)
ylabel("test accuracy")
title("held-out accuracy, same split")

knn: 0.983
tree: 0.95
forest: 0.967
bayes: 0.917
· bar: figure 10 (4 points)

Accuracy is computed here rather than taken from .score because naive_bayes defines predict and not score, and a comparison where one model is measured differently from the others is not a comparison. The boundary is a straight line, which is why the two that can draw one do best and the forest's extra capacity buys nothing.

The remaining constructors round out the table: two more clustering handles, a mixture model, a boosted regressor and two more binary classifiers.

km = kmedians_model(G, 3, seed = 5)
kd = kmedoids_model(G, 3, seed = 5)
gm = gmm_model(G, 3, seed = 5)
print("kmedians label for row 0: {km.predict(G)[0]}")
print("kmedoids label for row 0: {kd.predict(G)[0]}")
print("gmm log-likelihood: {round(gm.log_likelihood, 1)}")

gb = gradient_boosting_model(X, y, 20, learning_rate = 0.1, max_depth = 3)
print("gradient boosting r2: {round(gb.score(X, y), 3)}")

sv = svm_model(Xc, yc, kernel = "linear", C = 1.0)
lg = logistic_model(Xc, yc)
print("svm accuracy:      {round(sv.score(Xc, yc), 3)}")
print("logistic accuracy: {round(lg.score(Xc, yc), 3)}")

svr = svr_model(X, y)
print("svr r2: {round(svr.score(X, y), 3)}")

kmedians label for row 0: 1
kmedoids label for row 0: 1
gmm log-likelihood: [-448.7, -448.7, -448.7, -448.7]
gradient boosting r2: 0.953
svm accuracy:      0.995
logistic accuracy: 0.995
svr r2: 0.596

kmedoids_model picks its centres from among the rows it was given, so kd.predict can only ever return a label matching an actual data point's neighbourhood -- never an averaged one the way kmeans_model would.

Overloads: tree_model

Case: regression

kind = "regression" splits each node by variance reduction and predicts the mean of a leaf's training targets -- a continuous number.

Xtr = randn(80, 2, seed = 53)
ytr = Xtr[:, 0] * 3 + Xtr[:, 1] * -2 + 0.1 * randn(80, seed = 55)
Xte = randn(20, 2, seed = 54)
yte = Xte[:, 0] * 3 + Xte[:, 1] * -2
trm = tree_model(Xtr, ytr, max_depth = 4, kind = "regression")
print("tree_model regression R2 on new data: {round(trm.score(Xte, yte), 3)}")

tree_model regression R2 on new data: 0.799

Case: classification

kind = "classification" splits each node by Gini impurity instead and predicts the majority class among a leaf's training labels.

Xtrc = randn(80, 2, seed = 56)
ytrc = zeros(80)
for i in 0 to 79
    if Xtrc[i, 0] + Xtrc[i, 1] > 0
        ytrc[i] = 1
    end if
end for
Xtec = randn(20, 2, seed = 57)
ytec = zeros(20)
for i in 0 to 19
    if Xtec[i, 0] + Xtec[i, 1] > 0
        ytec[i] = 1
    end if
end for
tcm = tree_model(Xtrc, ytrc, max_depth = 4, kind = "classification")
print("tree_model classification accuracy on new data: {round(tcm.score(Xtec, ytec), 3)}")

tree_model classification accuracy on new data: 0.9

random_forest_model takes the same kind= and branches the same way, tree by tree.

Overloads: svm_model

Case: linear

kernel = "linear" fits a straight decision boundary -- fast, but blind to a relationship that isn't linearly separable.

Xx = [0.0, 0.0; 0.0, 1.0; 1.0, 0.0; 1.0, 1.0; 0.1, 0.1; 0.1, 0.9; 0.9, 0.1; 0.9, 0.9]
Yx = [0.0, 1.0, 1.0, 0.0, 0.0, 1.0, 1.0, 0.0]
svl = svm_model(Xx, Yx, kernel = "linear", C = 1.0)
print("linear kernel accuracy on XOR: {round(svl.score(Xx, Yx), 3)}")

linear kernel accuracy on XOR: 0.5

Case: rbf

kernel = "rbf" maps into an implicit, much higher-dimensional space through a Gaussian kernel scaled by gamma, so it can separate a boundary no straight line can.

Xx = [0.0, 0.0; 0.0, 1.0; 1.0, 0.0; 1.0, 1.0; 0.1, 0.1; 0.1, 0.9; 0.9, 0.1; 0.9, 0.9]
Yx = [0.0, 1.0, 1.0, 0.0, 0.0, 1.0, 1.0, 0.0]
svr = svm_model(Xx, Yx, kernel = "rbf", gamma = 2.0, C = 1.0)
print("rbf kernel accuracy on XOR:    {round(svr.score(Xx, Yx), 3)}")

rbf kernel accuracy on XOR:    1

Case: poly

kernel = "poly" raises the dot product to degree (shifted by coef0) -- a middle ground between "linear" and "rbf".

Xx = [0.0, 0.0; 0.0, 1.0; 1.0, 0.0; 1.0, 1.0; 0.1, 0.1; 0.1, 0.9; 0.9, 0.1; 0.9, 0.9]
Yx = [0.0, 1.0, 1.0, 0.0, 0.0, 1.0, 1.0, 0.0]
svp = svm_model(Xx, Yx, kernel = "poly", degree = 3, coef0 = 1.0, C = 1.0)
print("poly kernel accuracy on XOR:   {round(svp.score(Xx, Yx), 3)}")

poly kernel accuracy on XOR:   0.75

Overloads: knn_model

Case: regression

kind = "regression" predicts the mean of the k nearest neighbours' target values -- a number that can fall between any two training targets.

Xk = [0.0; 1.0; 2.0; 3.0]
yk = [10.0, 12.0, 30.0, 32.0]
knr = knn_model(Xk, yk, 2, kind = "regression")
q = [1.5] as matrix(1, 1)
print("knn_model regression predict(1.5): {knr.predict(q)}")   # mean of the 2 nearest y's

knn_model regression predict(1.5): [21]

Case: classification

kind = "classification" predicts the majority label among the k nearest neighbours instead -- always one of the training labels, never an average.

Xkc = [0.0; 1.0; 2.0; 3.0]
ykc = [0, 0, 1, 1]
knc = knn_model(Xkc, ykc, 3, kind = "classification")
q = [1.5] as matrix(1, 1)
print("knn_model classification predict(1.5): {knc.predict(q)}")   # majority vote among 3 nearest labels

knn_model classification predict(1.5): [0]

The shared verbs

score needs y for a supervised kind and refuses it for a clustering one; the GMM returns a total log-likelihood, not scikit-learn's mean-per-sample. Stages are named by string, the same convention pmap and spawn use.

Overloads: predict

Case: supervised or clustering model

model.predict(Xnew) returns a length-M vector -- one prediction per row of Xnew.

Xp = randn(30, 2, seed = 71)
yp = Xp * [1.0, -1.0] + 0.05 * randn(30, seed = 72)
mp = ols_model(Xp, yp)
print("model .predict returns a prediction vector: {round(mp.predict(Xp[0:1, :]), 3)}")

model .predict returns a prediction vector: [-1.1, -1.639]

Case: estimation filter state

state.predict(...) on a Kalman/EKF/UKF/particle state instead returns a whole new state Record -- the time-update step, not a prediction vector. The Estimation section below covers the filter-specific argument shapes.

s0 = kalman_init([0.0, 0.0], [1.0, 0, 0, 1.0] as matrix(2, 2))
F = [1, 0, 1, 1] as matrix(2, 2)
Q = [0.01, 0, 0, 0.01] as matrix(2, 2)
s1 = s0.predict(F, Q)
print("filter .predict returns a new state Record; its x field: {round(s1.x, 3)}")

filter .predict returns a new state Record; its x field: [0, 0]

Overloads: score

Case: regression

On a regression model (ols_model, ridge_model, svr_model, ...), .score(X, y) returns R².

Xr = randn(60, 2, seed = 61)
beta = [1.5, -0.7]
yr = Xr * beta + 0.05 * randn(60, seed = 62)
om = ols_model(Xr, yr)
print("regression score is R2: {round(om.score(Xr, yr), 4)}")

regression score is R2: 0.999

Case: classification

On a classifier (logistic_model, tree_model(kind="classification"), ...), .score(X, y) returns accuracy instead.

Xc2 = randn(60, 2, seed = 63)
yc2 = zeros(60)
for i in 0 to 59
    if Xc2[i, 0] + Xc2[i, 1] > 0
        yc2[i] = 1
    end if
end for
lgm = logistic_model(Xc2, yc2)
print("classification score is accuracy: {round(lgm.score(Xc2, yc2), 4)}")

classification score is accuracy: 1

Case: k-means

On a kmeans_model, .score(X) takes no y and returns negative inertia (the negated within-cluster sum of squared distances -- negative so that "higher is better" still holds).

Gs = vstack(randn(20, 2, seed = 64) + 4, randn(20, 2, seed = 65) - 4)
km = kmeans_model(Gs, 2, seed = 5)
print("kmeans score is negative inertia: {round(km.score(Gs), 4)}")

kmeans score is negative inertia: -59.3334

Case: PCA

On a pca_model, .score(X) returns the summed explained-variance fraction of the fitted components.

Gs = vstack(randn(20, 2, seed = 64) + 4, randn(20, 2, seed = 65) - 4)
pm = pca_model(Gs, 1)
print("pca score is summed explained variance: {round(pm.score(Gs), 4)}")

pca score is summed explained variance: 0.9807

Case: GMM

On a gmm_model, .score(X) returns the total log-likelihood of X under the fitted mixture.

Gs = vstack(randn(20, 2, seed = 64) + 4, randn(20, 2, seed = 65) - 4)
gm = gmm_model(Gs, 2, seed = 5)
print("gmm score is total log-likelihood: {round(gm.score(Gs), 2)}")

gmm score is total log-likelihood: -127.81

Recurrent cells: LSTM

Part of the tracked-tensor family with dense, relu, conv2d and gru_cell: a weight passed through param(...) gets a real gradient from grad(loss, wrt=...), backpropagated through time by the ordinary tape. No separate recurrent machinery.

lstm_forward is what you reach for when you have a whole sequence up front and just want the hidden states out the other end; it is also how gradients flow back through time here — because every step reuses the *same* tracked weight Tensors, grad sees the whole unrolled computation as one graph and backpropagates through every timestep in one call. The example below runs a 3-step sequence through a freshly initialized LSTM, computes a toy loss on the last hidden state, and takes the gradient with respect to one of the twelve weight matrices:

raw = lstm_init(2, 3, seed = 1)
w = {wf=param(raw.wf), uf=param(raw.uf), bf=param(raw.bf),
     wi=param(raw.wi), ui=param(raw.ui), bi=param(raw.bi),
     wo=param(raw.wo), uo=param(raw.uo), bo=param(raw.bo),
     wg=param(raw.wg), ug=param(raw.ug), bg=param(raw.bg)}
Xseq = [0.5, -0.3, 0.8; 0.2, 0.4, -0.1]   # (input_size, seq_len)

hs = lstm_forward(Xseq, [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], w)
h_last = hs[length(hs) - 1]
loss = sum(h_last .^ 2)
g = grad(loss, wrt = w.wf)                # through time, on the same tape

print("loss: {round(loss, 4)}")
print("dL/dwf: {round(g, 4)}")

loss: 0.0027
dL/dwf: [0.0001, 0; 0.0005, 0; 0.0021, -0.0002]

lstm_cell is the single timestep lstm_forward unrolls — useful on its own when you are stepping a sequence by hand, e.g. inside a control loop that needs to look at h/c between steps:

w1 = {wf=0.5, uf=0.5, bf=0.0, wi=0.3, ui=0.2, bi=0.0,
      wo=0.4, uo=0.1, bo=0.0, wg=0.6, ug=0.3, bg=0.0}
hc = lstm_cell(1.0, 0.0, 0.0, w1)
print("h = {round(hc.h, 4)}, c = {round(hc.c, 4)}")

h = 0.1791, c = 0.3085

Recurrent cells: GRU

The gated recurrent unit: fewer gates than an LSTM (one hidden state, not two), so gru_cell returns h directly rather than a {h, c} record.

The example below mirrors the LSTM one: unroll a GRU over a short sequence with gru_forward, take a gradient of a toy loss on the last hidden state back through every timestep, then reset the tape and show gru_cell stepping a single timestep by hand:

raw = gru_init(2, 3, seed = 1)
w = {wz=param(raw.wz), uz=param(raw.uz), bz=param(raw.bz),
     wr=param(raw.wr), ur=param(raw.ur), br=param(raw.br),
     wh=param(raw.wh), uh=param(raw.uh), bh=param(raw.bh)}
Xseq = [0.5, -0.3, 0.8; 0.2, 0.4, -0.1]   # (input_size, seq_len)

hs = gru_forward(Xseq, [0.0, 0.0, 0.0], w)
h_last = hs[length(hs) - 1]
print("gru_forward last hidden state: {round(h_last, 4)}")

loss = sum(h_last .^ 2)
g = grad(loss, wrt = w.wz)
print("dL/dwz: {round(g, 4)}")
tape_reset()

h1 = gru_cell([1.0], [1.0], {wz=0.5, uz=0.5, bz=0.0, wr=0.3, ur=0.3, br=0.0, wh=0.2, uh=0.2, bh=0.0})
print("gru_cell one step: {round(h1, 4)}")

gru_forward last hidden state: [-0.542; -0.0036; -0.0815]
dL/dwz: [0.2321, -0.0486; -0.0006, 0.0001; -0.0079, 0.0034]
gru_cell one step: [0.5012]

Local LLM inference (experimental)

Only in builds with --features llm. Runs on-device, CPU only, no API key. This build of the engine was compiled without that feature, so the two functions below are documented but not exercised by a runnable example here — calling either one on a plain build fails with a message naming the missing feature flag, the same pattern save_model/load_model use for --features h5-models further down this chapter.

Feedforward layers and activations

The plain building blocks every network above is made from. All of them are tracked-tensor aware: pass a param(...)-wrapped weight and grad sees through them.

The next example threads a small computation through every piece here at once: a dense layer feeds relu then sigmoid, grad differentiates a toy loss with respect to the weight matrix, tape_reset clears that computation before an unrelated softmax_rows call, and a final block shows dropout zeroing roughly half of eight ones and rescaling the rest.

w = param([0.5, -0.3; 0.2, 0.8])
b = param([0.1, -0.1])
x = [1.0, 2.0; -0.5, 1.5]      # two examples, two features each

z = dense(x, w, b)
a = relu(z)
p = sigmoid(a)
loss = sum(p .^ 2)
g = grad(loss, wrt = w)
print("relu:    {round(a, 4)}")
print("sigmoid: {round(p, 4)}")
print("dL/dw:   {round(g, 4)}")
tape_reset()                    # clear the tape before the next computation


scores = [2.0, 1.0, 0.1; 0.5, 0.5, 0.5]
probs = softmax_rows(scores)
print("softmax row 0: {round(probs[0, :], 4)}")

seed(3)
kept = dropout(ones(1, 8), 0.5)
print("dropout output: {kept}")
print("survivors:      {sum(kept > 0)} of 8, each scaled by 1/(1-rate)")

relu:    [1, 1.4; 0, 1.25]
sigmoid: [0.7311, 0.8022; 0.5, 0.7773]
dL/dw:   [0.2875, 0.12; 0.5749, 0.9128]
softmax row 0: [0.659, 0.2424, 0.0986]
dropout output: [2, 2, 0, 0, 2, 0, 2, 0]
survivors:      4 of 8, each scaled by 1/(1-rate)

Convolution and pooling

conv1d slides one learnable kernel over a signal; conv2d does the same over an image. maxpool2d/avgpool2d downsample by window.

The four functions below cover both dimensionalities you meet in practice: conv1d for a signal, conv2d for an image, and the two pooling functions for shrinking a feature map afterward. The example convolves a short signal with an edge-detecting kernel, convolves a small image with a diagonal-edge kernel, and then shows maxpool2d/avgpool2d shrinking the same image by taking the largest and the average value in each 2x2 block.

x1d = [1.0, 2.0, 3.0, 4.0, 5.0]
edge_kernel = [1.0, 0.0, -1.0]
y1d = conv1d(x1d, edge_kernel)
print("conv1d edge response: {y1d}")

img = [1, 5, 2, 8; 4, 3, 9, 1; 7, 2, 6, 0; 1, 1, 1, 1]
k2d = [1, 0; 0, -1]
y2d = conv2d(img, k2d, 1, "valid")
print("conv2d output:\n{y2d}")

print("maxpool2d: {maxpool2d(img, 2)}")
print("avgpool2d: {avgpool2d(img, 2)}")

conv1d edge response: [-2, -2, -2]
conv2d output:
[-2, -4, 1; 2, -3, 9; 6, 1, 5]
maxpool2d: [5, 9; 7, 6]
avgpool2d: [3.25, 5; 2.75, 2]

Layer normalization and attention

The transformer building blocks: layer_norm normalizes each row of its input independently (unlike batch norm, it needs no other rows in the batch to work); scaled_dot_product_attention is the core attention computation; multi_head_attention runs several of those in parallel over learned projections; transformer_block wraps attention and a feed-forward layer, each behind a residual connection and a layer norm.

The example below builds up from the smallest piece to the whole block: layer_norm on a two-row matrix with very different scales per row shows that each row is normalized independently of the other; a hand-built q/k/v triple runs through scaled_dot_product_attention directly; wrapping the identity matrix as a single attention head shows that multi_head_attention with one identity head reduces to plain self-attention, which is a useful sanity check when you are debugging a custom set of head weights; and a full transformer_block runs the same q through attention, a feed-forward layer, and two layer norms in one call.

x = [1.0, 2.0, 3.0; 10.0, 20.0, 30.0]
normed = layer_norm(x, [1.0, 1.0, 1.0], [0.0, 0.0, 0.0])
print("layer_norm (scale-invariant per row): {normed}")

q = [1.0, 0.5; -0.5, 1.0; 0.2, -0.3]
k = [1.0, 0.0; 0.0, 1.0; 0.5, 0.5]
v = [1.0, 2.0; 3.0, 4.0; 5.0, 6.0]
attn = scaled_dot_product_attention(q, k, v)
print("attention output: {round(attn, 3)}")

eye2 = [1.0, 0.0; 0.0, 1.0]
head = {wq = eye2, wk = eye2, wv = eye2, wo = eye2}
mh = multi_head_attention(q, (head,), 1)
plain = scaled_dot_product_attention(q, q, q)
print("multi_head_attention with 1 identity head matches plain self-attention: {sum(abs(mh - plain)) < 1e-9}")

weights = {heads = (head,), ln1_gamma = [1.0, 1.0], ln1_beta = [0.0, 0.0],
           w1 = [0.2, -0.1; 0.1, 0.3], b1 = [0.0, 0.0],
           w2 = [0.1, -0.2; 0.3, 0.1], b2 = [0.0, 0.0],
           ln2_gamma = [1.0, 1.0], ln2_beta = [0.0, 0.0]}
block_out = transformer_block(q, weights)
print("transformer_block: {round(block_out, 3)}")

layer_norm (scale-invariant per row): [-1.224736, 0, 1.224736; -1.224745, 0, 1.224745]
attention output: [2.872, 3.872; 3.25, 4.25; 2.872, 3.872]
multi_head_attention with 1 identity head matches plain self-attention: true
transformer_block: [1, -1; -1, 1; 1, -1]

Training loops and optimizers

Three ways to train, from least to most control. quick_mlp fits a whole network in one call. compile/.fit builds a sequential pipeline and trains it with a chosen optimizer. train_loop/optimizer_step run underneath both, and are also callable directly when you are writing the loop yourself.

sgd(w, 0.01)/adam(w, lr=0.01) take the tensor to optimize as their first argument; called with no positional argument at all (adam(lr=0.01)) they instead return a spec-only "recipe" for optimizer=. Every other update rule below plugs into the exact same optimizer_step(opt, grad) call — only the constructor changes:

x = [1.0, 2.0, 3.0, 4.0]
y_true = [2.0, 4.0, 6.0, 8.0]
w = param(0.0)
function loss_fn(w)
    pred = w * x
    diff = pred - y_true
    return mean(diff .^ 2)
end function
opt = sgd(w, 0.01)
losses = train_loop("loss_fn", w, 100, optimizer = "opt")
print("train_loop (sgd): loss {round(losses[0], 3)} -> {round(losses[99], 5)}")
tape_reset()

w2v = 0.0
opt2 = adam(param(w2v), lr = 0.1)
for t = 0 to 19
    tape_reset()
    w2 = param(w2v)
    pred = w2 * x
    loss2 = mean((pred - y_true) .^ 2)
    g = grad(loss2, wrt = w2)
    opt2 = optimizer_step(opt2, g)
    w2v = opt2._params[0]
end for
print("optimizer_step (adam): w -> {round(w2v, 4)} (target 2.0)")
tape_reset()

# One step each from w=1.0, grad=3.0 -- same optimizer_step(opt, grad) call.

w3 = param(1.0)
grad3 = [3.0]
r1 = optimizer_step(nesterov_sgd(w3, lr = 0.05, momentum = 0.9), grad3)
r2 = optimizer_step(adamw(w3, lr = 0.05), grad3)
r3 = optimizer_step(nadam(w3, lr = 0.05), grad3)
r4 = optimizer_step(adamax(w3, lr = 0.05), grad3)
r5 = optimizer_step(rmsprop(w3, lr = 0.05), grad3)
r6 = optimizer_step(adagrad(w3, lr = 0.05), grad3)
r7 = optimizer_step(adadelta(w3), grad3)
print("nesterov_sgd -> {round(r1._params[0], 4)}, adamw -> {round(r2._params[0], 4)}, nadam -> {round(r3._params[0], 4)}")
print("adamax -> {round(r4._params[0], 4)}, rmsprop -> {round(r5._params[0], 4)}, adagrad -> {round(r6._params[0], 4)}, adadelta -> {round(r7._params[0], 4)}")

Xxor = [0.0, 0.0; 0.0, 1.0; 1.0, 0.0; 1.0, 1.0]
Yxor = [0.0; 1.0; 1.0; 2.0]
net = input(2) |> dense(4, activation = "relu") |> dense(1, activation = "none") |>
      compile(optimizer = lr_plateau(lr = 0.1, patience = 5, base = "adam"), loss = "mse", seed = 7)
net = net.fit(Xxor, Yxor, epochs = 60)
n = length(net.loss_history)
print("compile + fit(lr_plateau): loss {round(net.loss_history[0], 4)} -> {round(net.loss_history[n - 1], 4)}")

net2 = input(2) |> dense(4, activation = "relu") |> dense(1, activation = "none") |>
       compile(optimizer = lr_adaptive(lr = 0.05), loss = "mse", seed = 7)
net2 = net2.fit(Xxor, Yxor, epochs = 60)
n2 = length(net2.loss_history)
print("compile + fit(lr_adaptive): loss {round(net2.loss_history[0], 4)} -> {round(net2.loss_history[n2 - 1], 4)}")

quick = quick_mlp(Xxor, Yxor, seed = 3)
nq = length(quick.loss_history)
print("quick_mlp: loss {round(quick.loss_history[0], 4)} -> {round(quick.loss_history[nq - 1], 4)}")

train_loop (sgd): loss 30 -> 0
optimizer_step (adam): w -> 1.7775 (target 2.0)
nesterov_sgd -> 0.715, adamw -> 0.9495, nadam -> 0.9263
adamax -> 0.95, rmsprop -> 0.8419, adagrad -> 0.95, adadelta -> 0.9996
compile + fit(lr_plateau): loss 3.3387 -> 0.0341
compile + fit(lr_adaptive): loss 0.3897 -> 0.0002
quick_mlp: loss 3.1877 -> 0.0922

Overloads: sgd

nesterov_sgd, adam, adamw, nadam, adamax, rmsprop, adagrad and adadelta all share this same two-call-shape overload.

wv = param(1.0)
bound = sgd(wv, lr = 0.1)
print("bound optimizer's _params: {bound._params}")   # [1] -- ready for optimizer_step


recipe = sgd(lr = 0.1)
print(recipe)                                          # model(sgd) {lr} -- no _params yet

bound optimizer's _params: [1]
model(sgd) {lr = 0.1}

Building the network: layer specs and sequential

The input(2) |> dense(4, activation = "relu") |> ... |> compile(...) pipeline used above is sugar over a smaller set of primitives: dense_layer and dropout_layer build plain, unfitted layer-spec Records (no weights yet), and sequential takes a list of them and returns an initialized — but not yet trained — kind="sequential" Model, exactly the same shape compile itself returns. mlp_classifier is a shortcut over the same two steps for the common "stack of dense+relu layers ending in raw class logits" case.

The example below builds the *same kind* of network two different ways. First, directly from primitives: two dense_layer calls with a dropout_layer in between, assembled by sequential, then fitted like any other kind="sequential" model. Second, via mlp_classifier, which skips writing out the layer list by hand for the common classification shape — its output is raw per-class logits, so a real classifier would still run them through softmax_rows before reading off a class.

dl1 = dense_layer(2, 8, activation = "relu")
dl2 = dropout_layer(0.3)
dl3 = dense_layer(8, 1, activation = "none")
net = sequential((dl1, dl2, dl3), seed = 3)
print("layer count: {length(net.layers)}")

Xxor = [0.0, 0.0; 0.0, 1.0; 1.0, 0.0; 1.0, 1.0]
Yxor = [0.0; 1.0; 1.0; 2.0]
net = net.fit(Xxor, Yxor, epochs = 60, optimizer = adam(lr = 0.05), loss = "mse")
n = length(net.loss_history)
print("sequential + dropout_layer: loss {round(net.loss_history[0], 4)} -> {round(net.loss_history[n - 1], 4)}")

mc = mlp_classifier(2, [8], 2, seed = 3)
Yc = [0, 1, 1, 0]
mc = mc.fit(Xxor, Yc, epochs = 60, loss = "cross_entropy")
nc = length(mc.loss_history)
print("mlp_classifier: loss {round(mc.loss_history[0], 4)} -> {round(mc.loss_history[nc - 1], 4)}")
preds = mc.predict(Xxor)
print("mlp_classifier raw logits shape: {rows(preds)} x {cols(preds)}")

layer count: 3
sequential + dropout_layer: loss 3.0224 -> 0.5
mlp_classifier: loss 0.7146 -> 0.6938
mlp_classifier raw logits shape: 4 x 2

Ready-made architectures

One call each fits a whole network, for when the question is about your data rather than about the network.

The two examples below build the smallest possible dataset for each architecture — two hand-drawn 6x6 images with different dot patterns for simple_cnn, and two short rising/falling sequences for simple_rnn_classifier — fit for a handful of epochs, and check that .predict comes back with one row per input sample:

net = simple_cnn([6, 6], 2, seed = 3)
img1 = zeros(6, 6)
img1[2, 2] = 1.0
img1[2, 3] = 1.0
img1[3, 2] = 1.0
img1[3, 3] = 1.0
img2 = zeros(6, 6)
img2[1, 2] = 1.0
img2[2, 1] = 1.0
img2[2, 3] = 1.0
img2[3, 2] = 1.0
net = net.fit((img1, img2), [0, 1], epochs = 5, lr = 0.1)
print("simple_cnn predict rows: {rows(net.predict((img1, img2)))}")

rnn = simple_rnn_classifier(1, 6, 2, seed = 5)
seq1 = [0.1, 0.3, 0.5, 0.7, 0.9] as matrix(1, 5)
seq2 = [0.9, 0.7, 0.5, 0.3, 0.1] as matrix(1, 5)
rnn = rnn.fit((seq1, seq2), [0, 1], epochs = 5, lr = 0.1)
print("simple_rnn_classifier predict rows: {rows(rnn.predict((seq1, seq2)))}")

simple_cnn predict rows: 2
simple_rnn_classifier predict rows: 2

Reinforcement learning: gridworld

A square grid, four actions, tabular temporal-difference learning.

The example below trains both algorithms on the same 4x4 grid for 300 episodes, compares their reward trend from the first 20 episodes to the last 20 (the real check that either one learned anything, since a constant reward with no crash proves nothing), then runs the greedy policy q_learning produced from the start state to the goal and counts the steps it takes:

env = gridworld_env(4)
start = env.reset()
print("start state: {start}")

q = q_learning(env, 300, alpha = 0.3, gamma = 0.9, epsilon = 0.1, seed = 7)
s = sarsa(env, 300, alpha = 0.3, gamma = 0.9, epsilon = 0.1, seed = 7)
print("q_learning reward: early {round(mean(q.rewards[0:19]), 2)} -> late {round(mean(q.rewards[280:299]), 2)}")
print("sarsa reward:      early {round(mean(s.rewards[0:19]), 2)} -> late {round(mean(s.rewards[280:299]), 2)}")

state = env.reset()
steps = 0
done = false
while done == false and steps < 20
    action = q.predict(state)
    result = env.step(state, action)
    state = result.next_state
    steps = steps + 1
    done = result.done
end while
print("greedy policy reached the goal in {steps} steps")

start state: 0
q_learning reward: early -4 -> late 4.45
sarsa reward:      early -1.5 -> late 4.6
greedy policy reached the goal in 6 steps

Estimation: Kalman, EKF, UKF, particle filters

The whole family goes through three verbs — predict, update, estimate — dispatching on the state's .kind and, within "kalman", on whether the second argument is a matrix (linear) or a function name (nonlinear, method="ekf" by default or "ukf"). Each call returns a new state; nothing is mutated.

Two things that cost time if you meet them the hard way. as matrix(r, c) fills column by column, so [1, 0, 1, 1] as matrix(2, 2) is the matrix whose first row is 1 1 — write the transition down and then read it back before trusting it. And jac= with method="ukf" is a rejected error rather than a silently ignored argument, because a UKF has no Jacobian to take.

The example below is the linear case — the one every other variant in the table above specializes: a constant-velocity [position, velocity] state tracked from noisy position-only measurements. Each loop iteration calls predict then update and keeps the new state, exactly the "returns a new state, nothing is mutated" contract every estimation function in this section follows; rmse at the end compares the raw noisy measurements against the filtered estimate, so you can see the filter is actually doing something rather than just running:

truth = zeros(60)
meas = zeros(60)
for t = 0 to 59
    truth[t] = 0.5 * t
    meas[t] = truth[t] + randn() * 2
end for

F = [1, 0, 1, 1] as matrix(2, 2)     # column-major: rows are [1, 1] and [0, 1]

Q = [0.01, 0, 0, 0.01] as matrix(2, 2)
H = [1, 0] as matrix(1, 2)
R = [4] as matrix(1, 1)

s = kalman_init([0, 0], [10, 0, 0, 10] as matrix(2, 2))
est = zeros(60)
for t = 0 to 59
    s = s.predict(F, Q)
    s = s.update(H, [meas[t]], R)
    est[t] = s.estimate()[0]
end for

print("measurement error: {round(rmse(truth, meas), 3)}")
print("filtered error:    {round(rmse(truth, est), 3)}")
scatter(0 to 59, meas, color = "#94a3b8", label = "measured")
plot(0 to 59, truth, color = "#0f172a", label = "truth")
plot(0 to 59, est, color = "#e11d48", width = 2, label = "Kalman")
legend()
title("constant-velocity Kalman filter")

measurement error: 2.197
filtered error:    1.193
· scatter: figure 11 (60 points)
· plot: figure 12 (60 points)
· plot: figure 13 (60 points)

Markov chains and hidden Markov models

Their own verbs — simulate/stationary, forward/viterbi — because neither model has an X to predict from.

The example below covers both models. First, a two-state weather chain: stationary gives the long-run fraction of time in each state, and chain.transition ^ 3 shows the 3-step transition probabilities computed with plain matrix power, no dedicated builtin needed. Second, a classic two-state HMM (healthy/sick, observing normal/cold/dizzy symptoms): forward totals the probability of one specific three-day observation sequence over every possible hidden-health-state path, and viterbi picks out the single most likely path instead.

P = [0.9, 0.5, 0.1, 0.5] as matrix(2, 2)   # column-major: rows [0.9 0.1], [0.5 0.5]

chain = markov_chain(P, 0)
print("stationary: {round(chain.stationary(), 4)}")
print("3-step:     {round(chain.transition ^ 3, 4)}")

trans = [0.7, 0.4, 0.3, 0.6] as matrix(2, 2)
emit  = [0.5, 0.1, 0.4, 0.3, 0.1, 0.6] as matrix(2, 3)
h = hmm(trans, emit, [0.6, 0.4])
print("P(normal, cold, dizzy) = {round(h.forward([0, 1, 2]), 5)}")
print("most likely path       = {h.viterbi([0, 1, 2])}")

stationary: [0.8333, 0.1667]
3-step:     [0.844, 0.156; 0.78, 0.22]
P(normal, cold, dizzy) = 0.03628
most likely path       = [0, 0, 1]

Computational geometry

A cell that reaches infinity has no finite vertex list to give you, so cells is empty or partial there and bounded is 0. On a small point set that is most of them. The example below computes the Voronoi diagram of a perfect square's four corners — every corner cell is unbounded, because with only four points there is nothing to close any of them off — and prints the one interior vertex that does show up: the square's centre, shared by all four unbounded cells.

pts = [0, 2, 0, 2, 0, 0, 2, 2] as matrix(4, 2)   # the corners of a square

v = voronoi(pts)
print("bounded: {v.bounded}")      # every corner cell is unbounded

print("cell 0:  {v.cells[0]}")     # the one finite vertex: the centre

bounded: [0, 0, 0, 0]
cell 0:  [1, 1]

Feature scaling

fit_scaler with .transform/.inverse_transform is the fit-once path — fit on training data, apply unchanged to a test set. The four one-shot verbs fit and apply in the same call, with nothing kept. All of them work column-wise on a matrix, whole on a vector, and per numeric column on a table.

One outlier is enough to show why the choice matters. normalize and standardize are the two everyday defaults — min-max to a fixed range, or z-score to mean 0 and unit variance — but both let a single extreme value dominate the whole scale, which is exactly what the example below shows by comparing standardize against the outlier-resistant robust_scale on the same eight numbers:

raw = [2, 4, 4, 4, 5, 5, 7, 900]
print("standardize:  {round(standardize(raw), 2)}")
print("robust_scale: {round(robust_scale(raw), 2)}")

standardize:  [-0.36, -0.35, -0.35, -0.35, -0.35, -0.35, -0.35, 2.47]
robust_scale: [-1.67, -0.33, -0.33, -0.33, 0.33, 0.33, 1.67, 597]

standardize squeezes every real value into a band a hundredth of a unit wide, because the outlier owns the standard deviation. robust_scale leaves them spread out and lets the outlier be the outlier.

quantile_normalize throws away the values entirely and keeps only the rank, so the outlier stops mattering at all:

print("quantile_normalize: {round(quantile_normalize(raw), 3)}")

quantile_normalize: [0, 0.143, 0.143, 0.143, 0.571, 0.571, 0.857, 1]

The four tied 4s and 5s land on the same rank, which is why the result has runs of equal values rather than eight evenly spaced ones.

Overloads: standardize

Case: vector

A plain Vec is standardized as a whole -- one mean, one std, across every element.

sv = [2.0, 4.0, 6.0, 8.0]
print("standardize on a vector: {round(standardize(sv), 3)}")

standardize on a vector: [-1.162, -0.387, 0.387, 1.162]

Case: matrix

An N-row, D-column Mat is standardized column by column -- each column gets its own mean and std, independent of the others.

sm = [1.0, 100.0; 2.0, 200.0; 3.0, 300.0]
print("standardize on a matrix, per column: {round(standardize(sm), 3)}")

standardize on a matrix, per column: [-1, -1; 0, 0; 1, 1]

Case: table

A table is standardized per numeric column, exactly like a matrix -- a non-numeric column (a label, say) passes through untouched instead of erroring.

st = table(a = [1.0, 2.0, 3.0], b = [10.0, 20.0, 30.0], label = ["x", "y", "z"])
print("standardize on a table, per numeric column (label untouched):")
print(standardize(st))

standardize on a table, per numeric column (label untouched):
 a   b  label
--  --  -----
-1  -1  x    
 0   0  y    
 1   1  z    

normalize, robust_scale and quantile_normalize all branch the same three ways.

Splitting and evaluation