  We propose two methods for generating non-Gaussian maps with fixed power
spectrum and bispectrum. The first makes use of a recently proposed rigorous,
non-perturbative, Bayesian framework for generating non-Gaussian distributions.
The second uses a simple superposition of Gaussian distributions. The former is
best suited for generating mildly non-Gaussian maps, and we discuss in detail
the limitations of this method. The latter is better suited for the opposite
situation, i.e. generating strongly non-Gaussian maps. The ensembles produced
are isotropic and the power spectrum can be jointly fixed; however we cannot
set to zero all other higher order cumulants (an unavoidable mathematical
obstruction). We briefly quantify the leakage into higher order moments present
in our method. We finally present an implementation of our code within the
HEALPIX package
