  We present ${\tt nimbus}$ : a hierarchical Bayesian framework to infer the
intrinsic luminosity parameters of kilonovae (KNe) associated with
gravitational-wave (GW) events, based purely on non-detections. This framework
makes use of GW 3-D distance information and electromagnetic upper limits from
multiple surveys for multiple events, and self-consistently accounts for finite
sky-coverage and probability of astrophysical origin. The framework is agnostic
to the brightness evolution assumed and can account for multiple
electromagnetic passbands simultaneously. Our analyses highlight the importance
of accounting for model selection effects, especially in the context of
non-detections. We show our methodology using a simple, two-parameter linear
brightness model, taking the follow-up of GW190425 with the Zwicky Transient
Facility (ZTF) as a single-event test case for two different prior choices of
model parameters -- (i) uniform/uninformative priors and (ii) astrophysical
priors based on surrogate models of Monte Carlo radiative transfer simulations
of KNe. We present results under the assumption that the KN is within the
searched region to demonstrate functionality and the importance of prior
choice. Our results show consistency with ${\tt simsurvey}$ -- an astronomical
survey simulation tool used previously in the literature to constrain the
population of KNe. While our results based on uniform priors strongly constrain
the parameter space, those based on astrophysical priors are largely
uninformative, highlighting the need for deeper constraints. Future studies
with multiple events having electromagnetic follow-up from multiple surveys
should make it possible to constrain the KN population further.
