# Multinomial Multinomial logistic family for unordered categorical responses. Usage ``` python Multinomial(n_categories) ``` [Multinomial](Multinomial.md#whittaker.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](OrderedCategorical.md#whittaker.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](Binomial.md#whittaker.Binomial), and for ordinal outcomes use [OrderedCategorical](OrderedCategorical.md#whittaker.OrderedCategorical). ## Parameters `n_categories: int` Number of response categories `K` (must be `>= 2`). Responses passed to [GAM.fit()](GAM.md#whittaker.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](OrderedCategorical.md#whittaker.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: ``` python 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](#category_intercepts) | Fitted per-category intercepts `alpha_1, ..., alpha_{K-1}`. | | [category_loadings](#category_loadings) | Fitted per-category loading coefficients `beta_1, ..., beta_{K-1}`. | | [n_categories](#n_categories) | Number of response categories `K`. | | [scale_known](#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`. `n_categories: int` ------------------------------------------------------------------------ ### scale_known Whether the dispersion/scale parameter is fixed rather than estimated. `scale_known: bool` ## Methods | Name | Description | |----|----| | [deviance()](#deviance) | Total multinomial deviance, implementing the family-specific `deviance` for | | [initialize()](#initialize) | Initialize category parameters and the starting linear predictor. | | [irls_update()](#irls_update) | Compute the P-IRLS working response and weights, implementing the family-specific | | [link()](#link) | Identity link, implementing the family-specific `link` for [Multinomial](Multinomial.md#whittaker.Multinomial). | | [link_derivative()](#link_derivative) | Derivative of the identity link, implementing the family-specific version for | | [link_inverse()](#link_inverse) | Identity inverse link, implementing the family-specific `link_inverse` for | | [log_likelihood()](#log_likelihood) | Multinomial log-likelihood, implementing the family-specific version for [Multinomial](Multinomial.md#whittaker.Multinomial). | | [simulate()](#simulate) | Draw random responses from the fitted category probabilities, implementing the | | [unit_deviance()](#unit_deviance) | Per-observation deviance contributions, implementing the family-specific version for | | [variance()](#variance) | Constant variance function, implementing the family-specific `variance` for | ------------------------------------------------------------------------ ### deviance() Total multinomial deviance, implementing the family-specific `deviance` for Usage ``` python deviance( y, mu, *, weights=None, ) ``` [Multinomial](Multinomial.md#whittaker.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. Usage ``` python initialize(y) ``` Implements the family-specific `initialize` for [Multinomial](Multinomial.md#whittaker.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 Usage ``` python irls_update( y, mu, eta, ) ``` update for [Multinomial](Multinomial.md#whittaker.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](Multinomial.md#whittaker.Multinomial). Usage ``` python link(mu) ``` 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 Usage ``` python link_derivative(mu) ``` [Multinomial](Multinomial.md#whittaker.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 Usage ``` python link_inverse(eta) ``` [Multinomial](Multinomial.md#whittaker.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](Multinomial.md#whittaker.Multinomial). Usage ``` python 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 Usage ``` python simulate( mu, scale, rng, ) ``` family-specific `simulate` for [Multinomial](Multinomial.md#whittaker.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 Usage ``` python unit_deviance( y, mu, ) ``` [Multinomial](Multinomial.md#whittaker.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 Usage ``` python variance(mu) ``` [Multinomial](Multinomial.md#whittaker.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`.