MLE Prediction Interface

Latent indices, mean-scale predictions, and labels after the MLE overhaul

This vignette documents the prediction contract used by the nonlinear MLE-style estimators after the 0.7 overhaul.

The contract

For the estimators where prediction makes sense, the public surface is now intentionally layered:

predict_lin(...): the latent linear index;
predict(...): the mean-scale or model-implied prediction object;
predict_label(...): only for classifiers, where hard labels are genuinely distinct from probabilities.

That contract avoids overloading predict() to sometimes mean probabilities and sometimes mean labels.

import numpy as np
import pandas as pd
import crabbymetrics as cm

np.set_printoptions(precision=4, suppress=True)

Binary logit

rng = np.random.default_rng(20260525)
n = 1200
x = rng.normal(size=(n, 3))
eta_true = -0.15 + x @ np.array([0.7, -0.45, 0.25])
p_true = 1.0 / (1.0 + np.exp(-eta_true))
y = rng.binomial(1, p_true).astype(np.int32)

logit = cm.Logit(max_iterations=300, gradient_tolerance=1e-8)
logit.fit(x, y)
print({key: logit.summary()[key] for key in ["converged", "iterations", "termination_reason", "objective"]})

check = x[:8]
logit_table = pd.DataFrame(
    {
        "eta": logit.predict_lin(check),
        "prob": logit.predict(check),
        "label@0.50": logit.predict_label(check),
        "label@0.35": logit.predict_label(check, cutoff=0.35),
    }
)
logit_table

{'converged': True, 'iterations': 8, 'termination_reason': 'Solver converged', 'objective': 763.2052710587117}

	eta	prob	label@0.50	label@0.35
0	0.499436	0.622327	1	1
1	1.248252	0.776997	1	1
2	-0.369529	0.408655	0	1
3	0.503944	0.623386	1	1
4	0.061744	0.515431	1	1
5	-0.327158	0.418932	0	1
6	-0.088640	0.477855	0	1
7	-0.765050	0.317551	0	0

The identity is exact by construction:

np.allclose(logit_table["prob"], 1.0 / (1.0 + np.exp(-logit_table["eta"])))

True

Multinomial logit

For multinomial logit, predict_lin() returns one classwise score per row, predict() returns class probabilities, and predict_label() returns the argmax class.

rng = np.random.default_rng(20260526)
n = 1500
coef = np.array(
    [
        [0.8, -0.2, 0.15],
        [-0.1, 0.65, -0.35],
        [0.25, -0.4, 0.45],
    ]
)
intercept = np.array([0.25, -0.15, 0.05])

x = rng.normal(size=(n, coef.shape[1]))
logits = x @ coef.T + intercept
logits = logits - logits.max(axis=1, keepdims=True)
probs = np.exp(logits)
probs = probs / probs.sum(axis=1, keepdims=True)
y = np.array([rng.choice(coef.shape[0], p=probs[i]) for i in range(n)], dtype=np.int32)

mlogit = cm.MultinomialLogit(max_iterations=300, gradient_tolerance=1e-8)
mlogit.fit(x, y)
print({key: mlogit.summary()[key] for key in ["converged", "iterations", "termination_reason", "objective"]})

check = x[:6]
eta_hat = mlogit.predict_lin(check)
prob_hat = mlogit.predict(check)
label_hat = mlogit.predict_label(check)

prob_table = pd.DataFrame(prob_hat, columns=["class0", "class1", "class2"])
prob_table["label"] = label_hat
prob_table

{'converged': True, 'iterations': 15, 'termination_reason': 'Solver converged', 'objective': 1420.4803365186447}

	class0	class1	class2	label
0	0.261540	0.584991	0.153469	1
1	0.316185	0.144446	0.539370	2
2	0.493774	0.365322	0.140903	0
3	0.369581	0.328764	0.301654	0
4	0.423588	0.164955	0.411456	0
5	0.238614	0.684788	0.076599	1

softmax = np.exp(eta_hat - eta_hat.max(axis=1, keepdims=True))
softmax = softmax / softmax.sum(axis=1, keepdims=True)
np.allclose(prob_hat, softmax)

True

Poisson

For Poisson, the latent index is the log conditional mean, and predict() exponentiates it back to the count scale.

rng = np.random.default_rng(20260527)
n = 900
x = rng.normal(size=(n, 2))
mu_true = np.exp(0.2 + x @ np.array([0.35, -0.3]))
y = rng.poisson(mu_true).astype(float)

poisson = cm.Poisson(max_iterations=250, tolerance=1e-8)
poisson.fit(x, y)
print({key: poisson.summary()[key] for key in ["converged", "iterations", "termination_reason", "objective"]})

check = x[:8]
poisson_table = pd.DataFrame(
    {
        "log_mu": poisson.predict_lin(check),
        "mu": poisson.predict(check),
    }
)
poisson_table

{'converged': True, 'iterations': 9, 'termination_reason': 'Solver converged', 'objective': 734.5201974599624}

	log_mu	mu
0	0.400609	1.492733
1	0.031761	1.032271
2	-0.337352	0.713657
3	-0.420010	0.657040
4	0.541645	1.718832
5	0.247413	1.280707
6	0.143784	1.154635
7	0.530882	1.700431

np.allclose(poisson_table["mu"], np.exp(poisson_table["log_mu"]))

True

Why this contract helps

The new split keeps three conceptually different objects separate:

the estimator’s linear score or index,
the model-implied prediction on the natural mean scale,
any final decision rule that thresholds or discretizes the mean-scale object.

That makes calibration work, diagnostics, and downstream scoring much cleaner than a classifier API where predict() skips straight to labels.