Multivariate Wallenius Distribution

• The multivariate Wallenius distribution is a generalized urn model that introduces category-specific weights to create biased sampling based on selection preferences.
• Likelihood-based and Bayesian methods, including ABC and the sphere-walk Metropolis algorithm, offer robust parameter estimation despite numerical challenges.
• Applications in preference analysis and contingency table modeling demonstrate the practical utility of the model in ranking categories and quantifying risk.

The multivariate Wallenius distribution, also known as the multivariate Wallenius noncentral hypergeometric distribution, generalizes the classical multivariate hypergeometric law by introducing positive, category-specific weights to the urn model of sampling without replacement. This framework is a principal tool for modelling preference or selection phenomena where categories (colours/groups) compete for representation and the sampling is biased according to a set of noncentrality parameters. Its formulation is analytically intractable for most practical problems but supports rigorous likelihood-based and Bayesian inference through carefully constructed computational procedures, including ABC and a novel sphere-walk Metropolis sampler (Grazian et al., 2017, Haines, 11 Nov 2025).

1. Urn Model and Probability Law

Consider an urn with $m$ colours, each corresponding to a category, where $N_i$ balls of colour $i$ (for $i=1,\dots,m$, $N_i > 0$) are present. $n$ balls are drawn sequentially, without replacement, with the probability of selecting colour $i$ at each draw proportional to its remaining count and a weight $w_i > 0$. These weights encode noncentrality or priority effects and are normalized so that $\sum_{i=1}^m w_i = 1$, placing $\mathbf{w}=(w_1,\ldots,w_m)$ in the interior of the $(m-1)$-simplex $\Delta^{m-1}$.

The joint probability mass function for a draw outcome $\mathbf{x}=(x_1,\dots,x_m)$ (with $\sum_i x_i = n$ and $0 \le x_i \le N_i$) is given by an integral form: $$f(\mathbf{x}\mid\mathbf{N},\mathbf{w}) = \left(\prod_{i=1}^m \binom{N_i}{x_i}\right)\frac{\int_0^1 \left[\prod_{i=1}^m(1-t^{w_i})^{N_i-x_i}\right] t^{\sum_{i=1}^m w_i x_i-1}\,\mathrm{d}t}{\mathrm{B}\left(\sum_i w_i x_i,\sum_i w_i(N_i-x_i)\right)},$$ where $\mathrm{B}(\alpha,\beta)$ is the Beta function. For $w_1 = \ldots = w_m$, the distribution reduces to the classical multivariate hypergeometric model, undoing the category-specific bias.

2. Inferential Procedures

Likelihood-Based Estimation

The log-likelihood for an observed outcome $\mathbf{x}$ is: $$\ell(\mathbf{w}; \mathbf{x}) = \sum_{i=1}^m \log\binom{N_i}{x_i} + \log\!\int_{0}^{1} \left[\prod_{i=1}^m (1-t^{w_i})^{N_i-x_i}\right] t^{\sum_i w_i x_i -1}\,\mathrm{d}t - \log\mathrm{B}\left(\sum_i w_i x_i, \sum_i w_i(N_i-x_i)\right).$$ Direct numerical maximization is complicated by the simplex constraint on $\mathbf{w}$. This is resolved by the additive log-ratio transform: $$\theta_i = \log\frac{w_i}{w_m},\quad i=1,\dots,m-1,$$ with inversion

$$w_i = \frac{\exp(\theta_i)}{1 + \sum_{j=1}^{m-1}\exp(\theta_j)},\qquad w_m = \frac{1}{1+\sum_{j=1}^{m-1}\exp(\theta_j)}.$$

Optimization employs quasi-Newton algorithms, using analytic or automatic differentiation. Under- and overflow in the numerical integration are controlled by extracting a pointwise exponent maximum $M(\mathbf{w})$ and scaling the integrand accordingly.

Bayesian Approaches

A uniform Dirichlet prior $\mathrm{Dir}(1,\dots,1)$ over $\Delta^{m-1}$ is standard. Posterior inference involves samples $\mathbf{w}^{(s)}$ with

$$\pi(\mathbf{w}|\mathbf{x}) \propto f(\mathbf{x}|\mathbf{w})\pi_0(\mathbf{w}).$$

Approximate Bayesian Computation (ABC) is employed when direct likelihood evaluation is prohibitive (Grazian et al., 2017). ABC rejection sampling proceeds by simulating preference profiles from candidate $\mathbf{w}'$, computing summary statistics (typically average relative frequencies per category)

$$\eta(\mathbf{x}) = \frac{1}{k} \sum_{h=1}^{k} (x_{h1}/n_h, \dots, x_{hc}/n_h),$$

and accepting when the variation-distance between simulated and observed summaries falls below a pre-specified tolerance $\varepsilon$.

The sphere-walk Metropolis (SWM) algorithm for MCMC sampling from the posterior (Haines, 11 Nov 2025) operates by embedding the simplex in Cartesian coordinates via an orthonormal basis, proposing moves along random directions on a sphere, and filtering in barycentric coordinates. Volume preservation (constant Jacobian) and symmetric proposal densities simplify the Metropolis acceptance ratio.

3. Numerical Challenges and Remedies

Likelihoods stemming from the Wallenius model can be extremely small; exponent scaling is critical for numerical stability. Integration over $[0,1]$ within the pmf is typically approximated using adaptive Gauss–Legendre or Gauss–Jacobi quadrature. In SWM, the choice of step-size $\delta$ is influential, with acceptance rates in the $20\%$–$40\%$ range optimal for $m$ up to $10$. Iterative re-projection onto $\Delta^{m-1}$ guards against drift outside the simplex. These procedures are essential to maintain computational feasibility in both frequentist and Bayesian analysis for moderate and large $m$.

4. Relation to Classical Models and Interpretability

The multivariate Wallenius distribution strictly generalizes the hypergeometric law. When all category weights are equal, $\mathbf{w}=(1/m,\dots,1/m)$, sampling is unbiased and reduces to the classical model. Increasing $\omega_i$ (or $w_i$) for a particular category favours its early selection. The weight vector thus encodes the interpretable rank or attractiveness of categories, foundational in applications to ranking, preference analysis, and contingency tables with fixed margins (Grazian et al., 2017, Haines, 11 Nov 2025).

5. Applications

Two substantive applications exemplify the model’s utility:

• Preference Data Modelling: In analysis of movie ratings, users label movies by genre and the Wallenius model is applied to infer genre popularity. ABC estimation of $\omega$ displays slight but statistically supported differences among genres, with Action and Sci-Fi genres evidencing higher preference weights. For Italian academic statisticians’ journal preferences, journals grouped into five categories yield a ranking with Methodology and Applied categories at the top, Probability at the bottom, and intermediate values for Computational and Econometrics. Pairwise probabilities $p_{ij} = \Pr(\omega_i > \omega_j)$ quantify uncertainties in these rankings (Grazian et al., 2017).
• Contingency Table Analysis: In bristle-fly chemical mortality data [Manly (1974)], each bristle-count group constitutes a colour, and the probability of death under chemical exposure is modelled by the Wallenius law. Maximum likelihood and Bayesian estimates of weights $w_i$ for each category quantify relative susceptibilities. Posterior summaries provide credible intervals for each $w_i$, offering interpretable measures of category-specific risk (Haines, 11 Nov 2025).

A summary table for bristle-fly data:

Bristle Count Group ($i$)	MLE $\widehat{w}_i$	Posterior Mean $w_i$
14	0.112	0.113
15	0.058	0.062
16	0.210	0.206
...	...	...
21–22	0.086	0.088

Interpretation: Larger $w_i$ for a group denotes higher susceptibility or “attractiveness” for selection.

6. Extensions and Prospects

The Wallenius model is currently focused on category-level analysis. An extension to nested hierarchical models permitting within-category heterogeneity is a plausible implication for future work. Efficient simulation routines exist (e.g., R package BiasedUrn). Open research directions include scalable computation for large $m,n$ and theoretical exploration of model identifiability and sensitivity to the choice of weight normalization and priors (Grazian et al., 2017, Haines, 11 Nov 2025).

7. Concluding Perspective

The multivariate Wallenius distribution enables principled inference on biased sampling schemes with analytic tractability deferred to numerical and simulation-based methods. Its capacity to generalize and interpret rank-driven selection phenomena, combined with robust computational algorithms for inference, establishes its relevance for advanced categorical data analysis, complex contingency tables, and structured preference data. Future developments in scalable hierarchical extensions and computation will further enhance its applicability to modern statistical modelling.