A realistic biomechanical thumb model would elucidate the functional consequences of orthopedic and neurological diseases and their treatments. We investigated whether a single parametric kinematic model can represent all thumbs, or whether different kinematic model structures are needed to represent different thumbs. We used Monte Carlo simulations to convert the anatomical variability in the kinematic model parameters into distributions of Denavit-Hartenberg parameters amenable for robotics-based models. Upon convergence (3550 simulations, where mean and coefficient of variance changed < 1% for the last 20 + % simulations) the distributions of Denavit-Hartenberg parameters appeared multimodal, in contrast to the reported unimodal distributions of the anatomy-based parameters. Cluster analysis and one-way analysis of variance confirmed four types of kinematic models (p < 0.0001) differentiated primarily by the biomechanically relevant order of MCP joint axes (in 65.2% of models, the flexion-extension axis was distal to the adduction-abduction axis); and secondarily by a detail specifying the direction of a common normal between successive axes of rotation. Importantly, this stochastic analysis of anatomical variability redefines the debate on whether a single generic biomechanical model can represent the entire population, or if subject-specific models are necessary. We suggest a practical third alternative: that anatomical and functional variability can be captured by a finite set of model-types.