This paper takes a Bayesian approach to the problem of dynamic system estimation, and illustrates how posterior densities for rather arbitrary system parameters or properties (such a frequency response, phase margin etc) may be numerically computed. In achieving this, the key idea of constructing an ergodic Markov chain with invariant distribution equal to the desired posterior is one borrowed from the mathematical statistics literature. An essential point of the work here is that, via the associated posterior computation from the Markov chain, error bounds on estimates are provided that do not rely on asymptotic in data length arguments, and hence they apply with arbitrary accuracy for arbitrarily short data records.