  We develop a framework for joint constraints on the CO luminosity function
based on power spectra (PS) and voxel intensity distributions (VID), and apply
this to simulations of COMAP, a CO intensity mapping experiment. This Bayesian
framework is based on a Markov chain Monte Carlo (MCMC) sampler coupled to a
Gaussian likelihood with a joint PS + VID covariance matrix computed from a
large number of fiducial simulations, and re-calibrated with a small number of
simulations per MCMC step. The simulations are based on dark matter halos from
fast peak patch simulations combined with the $L_\text{CO}(M_\text{halo})$
model of Li et al. (2016). We find that the relative power to constrain the CO
luminosity function depends on the luminosity range of interest. In particular,
the VID is more sensitive at both small and large luminosities, while the PS is
more sensitive at intermediate luminosities. The joint analysis is superior to
using either observable separately. When averaging over CO luminosities ranging
between $L_\text{CO} = 10^4-10^7L_\odot$, and over 10 cosmological realizations
of COMAP Phase 2, the uncertainties (in dex) are larger by 58 % and 30 % for
the PS and VID, respectively, when compared to the joint analysis (PS + VID).
This method is generally applicable to any other random field, with a
complicated likelihood, as long a fast simulation procedure is available.
