  We introduce a Bayesian solution to the problem of inferring the density
profile of strong gravitational lenses when the lens galaxy may contain
multiple dark or faint substructures. The source and lens models are based on a
superposition of an unknown number of non-negative basis functions (or "blobs")
whose form was chosen with speed as a primary criterion. The prior distribution
for the blobs' properties is specified hierarchically, so the mass function of
substructures is a natural output of the method. We use reversible jump Markov
Chain Monte Carlo (MCMC) within Diffusive Nested Sampling (DNS) to sample the
posterior distribution and evaluate the marginal likelihood of the model,
including the summation over the unknown number of blobs in the source and the
lens. We demonstrate the method on two simulated data sets: one with a single
substructure, and one with ten. We also apply the method to the g-band image of
the "Cosmic Horseshoe" system, and find evidence for more than zero
substructures. However, these have large spatial extent and probably only point
to misspecifications in the model (such as the shape of the smooth lens
component or the point spread function), which are difficult to guard against
in full generality.
