  We present the implementation of a Bayesian formalism within the
Stochastically Lighting Up Galaxies (SLUG) stellar population synthesis code,
which is designed to investigate variations in the initial mass function (IMF)
of star clusters. By comparing observed cluster photometry to large libraries
of clusters simulated with a continuously varying IMF, our formalism yields the
posterior probability distribution function (PDF) of the cluster mass, age, and
extinction, jointly with the parameters describing the IMF. We apply this
formalism to a sample of star clusters from the nearby galaxy NGC 628, for
which broad-band photometry in five filters is available as part of the Legacy
ExtraGalactic UV Survey (LEGUS). After allowing the upper-end slope of the IMF
($\alpha_3$) to vary, we recover PDFs for the mass, age, and extinction that
are broadly consistent with what is found when assuming an invariant Kroupa
IMF. However, the posterior PDF for $\alpha_3$ is very broad due to a strong
degeneracy with the cluster mass, and it is found to be sensitive to the choice
of priors, particularly on the cluster mass. We find only a modest improvement
in the constraining power of $\alpha_3$ when adding H$\alpha$ photometry from
the companion H$\alpha$-LEGUS survey. Conversely, H$\alpha$ photometry
significantly improves the age determination, reducing the frequency of
multi-modal PDFs. With the aid of mock clusters we quantify the degeneracy
between physical parameters, showing how constraints on the cluster mass that
are independent of photometry can be used to pin down the IMF properties of
star clusters.
