  We constrain the evolution of the rest-frame far-infrared (FIR) luminosity
function out to high redshift, by combining several pieces of complementary
information provided by the deep Balloon-borne Large-Aperture Submillimeter
Telescope surveys at 250, 350 and 500 micron, as well as other FIR and
millimetre data. Unlike most other phenomenological models, we characterise the
uncertainties in our fitted parameters using Monte Carlo Markov Chains. We use
a bivariate local luminosity function that depends only on FIR luminosity and
60-to-100 micron colour, along with a single library of galaxy spectral energy
distributions indexed by colour, and apply simple luminosity and density
evolution. We use the surface density of sources, Cosmic Infrared Background
(CIB) measurements and redshift distributions of bright sources, for which
identifications have been made, to constrain this model. The precise evolution
of the FIR luminosity function across this crucial range has eluded studies at
longer wavelengths (e.g., using SCUBA and MAMBO) and at shorter wavelengths
(e.g., Spitzer), and should provide a key piece of information required for the
study of galaxy evolution. Our adoption of Monte Carlo methods enables us not
only to find the best-fit evolution model, but also to explore correlations
between the fitted parameters. Our model-fitting approach allows us to focus on
sources of tension coming from the combination of data-sets. We specifically
find that our choice of parameterisation has difficulty fitting the combination
of CIB measurements and redshift distribution of sources near 1 mm. Existing
and future data sets will be able to dramatically improve the fits, as well as
break strong degeneracies among the models. [abridged]
