Inferring the number, rate and intrinsic properties of short gamma-ray bursts has been a long studied problem in the field. As it is closely related to the number of GW events expected for neutron star mergers, the topic has begun to be discussed int he literature again. However, the utilized techniques for GRBs still rely on improper statistical modeling V/Vmax estimators and in many cases, methods made for humor alone. I will discuss the use of Bayesian hierarchical models to infer population and object level parameters of inhomogeneous-Poisson process distributed populations. Techniques to handle high-dimensional selections effects will be introduced. The methodology will then be applied to sGRB population data with the aim of understand how many of these objects there are, where they are in the Universe and what are their properties under given modeling assumptions. The methodology is general, thus extensions to other populations can be made easily.