  Black hole (BH) superradiance can provide strong constraints on the
properties of ultralight bosons (ULBs). Since most of the previous work has
focused on the theoretical predictions, here we investigate the most suitable
statistical framework to constrain ULB masses and self-interactions. We argue
that a Bayesian approach provides a clear statistical interpretation, deals
with limitations regarding the reproducibility of existing BH analyses,
incorporates the full information from BH data, and allows us to include
additional nuisance parameters or to perform hierarchical modelling with BH
populations in the future. We demonstrate the feasibility of our approach using
mass and spin posterior samples for the X-ray binary BH M33 X-7 and, for the
first time in this context, the supermassive BH IRAS 09149-6206. We explain the
differences to existing ULB constraints in the literature and illustrate the
effects of various assumptions about the superradiance process (equilibrium
regime vs cloud collapse, higher occupation levels). As a result, our procedure
yields the most rigorous ULB constraints available in the literature, with
important implications for the QCD axion and axion-like particles. We encourage
all groups analysing BH data to publish likelihood functions or posterior
samples as supplementary material to facilitate this type of analysis.
