We propose a new method to infer the star formation histories of resolved stellar populations. With photometry one may plot observed stars on a colour-magnitude diagram (CMD) and then compare with synthetic CMDs representing different star formation histories. This has been accomplished hitherto by parametrising the model star formation history as a histogram, usually with the bin widths set by fixed increases in the logarithm of time. A best fit is then found with maximum likelihood methods and we consider the different means by which a likelihood can be calculated. We then apply Bayesian methods by parametrising the star formation history as an unknown number of Gaussian bursts with unknown parameters. This parametrisation automatically provides a smooth function of time. A Reversal Jump Markov Chain Monte Carlo method is then used to find both the most appropriate number of Gaussians, thus avoiding avoid overfitting, and the posterior probability distribution of the star formation rate. We apply our method to artificial populations and to observed data. We discuss the other advantages of the method: direct comparison of different parametrisations and the ability to calculate the probability that a given star is from a given Gaussian. This allows the investigation of possible sub-populations.