Most astronomers are implicit Bayesians.

An extension of the work of Loredo (1992) gives the joint probability distribution of model parameters for Poisson-distributed data in the presence of background.

I have an implementation in XSPEC that successfully finds the ``best-fit'' and allows prior distributions on parameters.

To get the distribution for one parameter we must integrate over the others. This is a difficult problem involving multi-dimensional integrals and I don't have a working solution at present (anyone want to help ?).