  Electromagnetic stimulation probes and modulates the neural systems that
control movement. Key to understanding their effects is the muscle recruitment
curve, which maps evoked potential size against stimulation intensity. Current
methods to estimate curve parameters require large samples; however, obtaining
these is often impractical due to experimental constraints. Here, we present a
hierarchical Bayesian framework that accounts for small samples, handles
outliers, simulates high-fidelity data, and returns a posterior distribution
over curve parameters that quantify estimation uncertainty. It uses a
rectified-logistic function that estimates motor threshold and outperforms
conventionally used sigmoidal alternatives in predictive performance, as
demonstrated through cross-validation. In simulations, our method outperforms
non-hierarchical models by reducing threshold estimation error on sparse data
and requires fewer participants to detect shifts in threshold compared to
frequentist testing. We present two common use cases involving electrical and
electromagnetic stimulation data and provide an open-source library for Python,
called hbMEP, for diverse applications.
