Multinomial logistic family for unordered categorical responses.
Multinomial(n_categories)
Multinomial models a categorical response with K unordered levels — for example, a choice among several unranked options — using a baseline-category logit model. Unlike OrderedCategorical, no assumption is made about the ordering of categories or a shared direction of covariate effects: each non-reference category k gets its own intercept alpha_k and its own loading coefficient beta_k that rescales the shared linear predictor eta, so different categories can respond differently (even in sign) to the same covariate effect. The final category K is fixed as the reference, with alpha_K = 0 and beta_K = 0. Use this family when the response is nominal (has no natural order) with more than two levels; for binary outcomes use Binomial, and for ordinal outcomes use OrderedCategorical.
Parameters
n_categories: int-
Number of response categories
K (must be >= 2). Responses passed to GAM.fit() should be integer-coded 1, 2, ..., K, with category K as the reference.
Notes
Category probabilities are obtained from a softmax over per-category logits built from the shared linear predictor eta:
P(Y = k \mid \eta) = \frac{\exp(\alpha_k + \beta_k \eta)}
{\sum_{j=1}^{K} \exp(\alpha_j + \beta_j \eta)},
\qquad \alpha_K = \beta_K = 0.
As with OrderedCategorical, this loss does not fit the standard GLM deviance framework: link and link_inverse are the identity on eta, and a custom irls_update drives P-IRLS while an inner maximum-likelihood step (_update_params) re-estimates the per-category intercepts alpha and loadings beta at each iteration. The reported deviance is -2 times the multinomial log-likelihood of the observed categories under the fitted probabilities,
D(y, \hat P) = -2 \sum_{i} \log \hat P(Y_i = y_i \mid \eta_i).
Examples
Fit a GAM to a three-level unordered categorical response:
import numpy as np
import whittaker as wk
rng = np.random.default_rng(0)
n = 300
x = np.linspace(-3, 3, n)
eta = np.sin(x)
alphas = np.array([0.0, 0.5])
betas = np.array([1.0, -1.5])
logits = np.column_stack(
[alphas[0] + betas[0] * eta, alphas[1] + betas[1] * eta, np.zeros(n)]
)
probs = np.exp(logits) / np.exp(logits).sum(axis=1, keepdims=True)
y = np.array([rng.choice([1, 2, 3], p=probs[i]) for i in range(n)], dtype=float)
data = {"x": x, "y": y}
model = wk.GAM("y ~ s(x)", family=wk.Multinomial(n_categories=3))
model.fit(data, method="REML")
print(model.summary())
GAM fit summary
============================================================
Formula: y ~ s(x)
Family: Multinomial(K=3)
Inference: REML
Observations: 300
Coefficients: 10
Parametric coefficients:
Term Estimate Std.Err z value p-value
------------------------ ---------- ---------- ---------- ----------
(Intercept) -0.0013 0.0231 -0.057 0.9545
Approximate significance of smooth terms:
Term EDF Ref.df Chi.sq p-value
------------------------ ------ ------ ---------- ----------
s(x) 5.39 6 87.175 1.169e-16
Total EDF: 6.39
Scale est: 1.000000
Deviance: 534.0030
Null dev: 4154.4017
Dev. expl: 87.1%
GCV score: 1.858347
AIC: 546.79
BIC: 570.46
Attributes
|
Name
|
Description
|
|
category_intercepts
|
Fitted per-category intercepts alpha_1, ..., alpha_{K-1}.
|
|
category_loadings
|
Fitted per-category loading coefficients beta_1, ..., beta_{K-1}.
|
|
n_categories
|
Number of response categories K.
|
|
scale_known
|
Whether the dispersion/scale parameter is fixed rather than estimated.
|
category_intercepts
Fitted per-category intercepts alpha_1, ..., alpha_{K-1}.
category_intercepts: NDArray | None
category_loadings
Fitted per-category loading coefficients beta_1, ..., beta_{K-1}.
category_loadings: NDArray | None
n_categories
Number of response categories K.
scale_known
Whether the dispersion/scale parameter is fixed rather than estimated.
Methods
|
Name
|
Description
|
|
deviance()
|
Total multinomial deviance, implementing the family-specific deviance for
|
|
initialize()
|
Initialize category parameters and the starting linear predictor.
|
|
irls_update()
|
Compute the P-IRLS working response and weights, implementing the family-specific
|
|
link()
|
Identity link, implementing the family-specific link for Multinomial.
|
|
link_derivative()
|
Derivative of the identity link, implementing the family-specific version for
|
|
link_inverse()
|
Identity inverse link, implementing the family-specific link_inverse for
|
|
log_likelihood()
|
Multinomial log-likelihood, implementing the family-specific version for Multinomial.
|
|
simulate()
|
Draw random responses from the fitted category probabilities, implementing the
|
|
unit_deviance()
|
Per-observation deviance contributions, implementing the family-specific version for
|
|
variance()
|
Constant variance function, implementing the family-specific variance for
|
deviance()
Total multinomial deviance, implementing the family-specific deviance for
deviance(
y,
mu,
*,
weights=None,
)
Multinomial.
Computes -2 times the multinomial log-likelihood of the observed categories y under the fitted category probabilities implied by the linear predictor eta = mu:
D(y, \hat P) = -2 \sum_{i=1}^{n} \log \hat P(Y_i = y_i \mid \eta_i).
Returns float(len(y)) as a placeholder if the family has not yet been fitted.
Parameters
y: NDArray-
Observed categories, coded 1, 2, ..., K.
mu: NDArray-
Linear predictor eta (this family uses the identity link, so mu and eta coincide).
weights: NDArray | None = None-
Accepted for interface compatibility; does not currently affect the computed deviance.
Returns
float-
The total deviance.
initialize()
Initialize category parameters and the starting linear predictor.
Implements the family-specific initialize for Multinomial: estimates starting intercepts alpha (from empirical log-odds relative to the reference category) and loadings beta (set to 1) via _init_params, then returns a starting linear predictor of all zeros for the outer P-IRLS loop to refine.
Parameters
y: NDArray-
Observed categories, coded
1, 2, ..., K, used to compute starting intercepts.
Returns
NDArray-
Initial linear predictor
eta, an array of zeros with the same shape as y.
irls_update()
Compute the P-IRLS working response and weights, implementing the family-specific
irls_update(
y,
mu,
eta,
)
update for Multinomial.
Overrides the default GLM IRLS step: because the multinomial deviance is not a standard exponential-family deviance in eta, this method first re-estimates the per-category intercepts alpha and loadings beta by maximum likelihood (via _update_params), initializing them on the first call (via _init_params). It then computes, for each observation i with observed category k = y_i - 1 and fitted probabilities p = P(Y_i = \cdot \mid \eta_i), the gradient and negative curvature of the multinomial log-likelihood with respect to eta_i:
\frac{\partial \ell_i}{\partial \eta_i} = \sum_{j=1}^{K-1} \beta_j \left(
\mathbb{1}[j = k] - p_j \right), \qquad
-\frac{\partial^2 \ell_i}{\partial \eta_i^2} = \sum_{j=1}^{K-1} \beta_j^2 p_j
(1 - p_j) - 2 \sum_{j < l} \beta_j \beta_l p_j p_l.
The working weight W is the (clipped) negative curvature, and the working response is the Newton step z = eta + grad / W; both feed into the outer P-IRLS smoothing loop.
Parameters
y: NDArray-
Observed categories, coded 1, 2, ..., K.
mu: NDArray-
Unused; present for interface compatibility with other families.
eta: NDArray-
Current linear predictor values.
Returns
tuple[NDArray, NDArray]-
The working response
z and working weights W, both of shape (n,).
link()
Identity link, implementing the family-specific link for Multinomial.
Because mu here already represents the shared linear predictor eta fed into the per-category softmax (rather than a mean response on the natural scale), the link is the identity: g(\mu) = \mu.
Parameters
mu: NDArray-
Linear predictor values.
Returns
NDArray-
The input
mu, unchanged.
link_derivative()
Derivative of the identity link, implementing the family-specific version for
Multinomial.
Since link is the identity, g'(\mu) = 1 everywhere.
Parameters
mu: NDArray-
Linear predictor values (only used to determine output shape).
Returns
NDArray-
Array of ones with the same shape as
mu.
link_inverse()
Identity inverse link, implementing the family-specific link_inverse for
Multinomial.
Complements link: since eta is passed straight through to _category_probs for the softmax computation, the inverse link is also the identity, g^{-1}(\eta) = \eta.
Parameters
eta: NDArray-
Linear predictor values.
Returns
NDArray-
The input
eta, unchanged.
log_likelihood()
Multinomial log-likelihood, implementing the family-specific version for Multinomial.
log_likelihood(
y,
mu,
scale,
*,
weights=None,
)
Recovers the log-likelihood from deviance via \ell = -D / 2, since deviance is defined as -2 times the log-likelihood.
Parameters
y: NDArray-
Observed categories, coded 1, 2, ..., K.
mu: NDArray-
Linear predictor eta.
scale: float-
Unused; present for interface compatibility with other families (the multinomial scale parameter is fixed, see scale_known).
weights: NDArray | None = None-
Forwarded to
deviance.
Returns
float-
The total log-likelihood of
y under the fitted model.
simulate()
Draw random responses from the fitted category probabilities, implementing the
simulate(
mu,
scale,
rng,
)
family-specific simulate for Multinomial.
For each observation, computes the category probabilities from the linear predictor eta = mu via _category_probs and draws one category from {1, ..., K} according to those probabilities using rng.choice.
Parameters
mu: NDArray-
Linear predictor eta at which to simulate.
scale: float-
Unused; present for interface compatibility with other families.
rng: object-
A random number generator exposing a
choice(a, p=...) method (e.g. a NumPy Generator).
Returns
NDArray-
Simulated categories, coded
1, 2, ..., K, one per row of mu.
Raises
RuntimeError-
If called before the family has been fitted.
unit_deviance()
Per-observation deviance contributions, implementing the family-specific version for
Multinomial.
For each observation i, returns -2 \log \hat P(Y_i = y_i \mid \eta_i), the per-observation contribution to deviance. Returns an array of ones if the family has not yet been fitted.
Parameters
y: NDArray-
Observed categories, coded 1, 2, ..., K.
mu: NDArray-
Linear predictor
eta (this family uses the identity link).
Returns
NDArray-
Array of shape
(n,) with the deviance contribution of each observation; these sum to the value returned by deviance.
variance()
Constant variance function, implementing the family-specific variance for
Multinomial.
The multinomial log-likelihood does not follow the mean-variance relationship used by standard GLM families; the working weights used by P-IRLS are instead derived directly from the curvature of the multinomial log-likelihood in irls_update. This method simply returns an array of ones so it is a no-op wherever a generic variance function might otherwise be referenced.
Parameters
mu: NDArray-
Linear predictor values (only used to determine output shape).
Returns
NDArray-
Array of ones with the same shape as
mu.