  We present an axisymmetric, equilibrium model for late-type galaxies which
consists of an exponential disk, a Sersic bulge, and a cuspy dark halo. The
model is specified by a phase space distribution function which, in turn,
depends on the integrals of motion. Bayesian statistics and the Markov Chain
Monte Carlo method are used to tailor the model to satisfy observational data
and theoretical constraints. By way of example, we construct a chain of 10^5
models for the Milky Way designed to fit a wide range of photometric and
kinematic observations. From this chain, we calculate the probability
distribution function of important Galactic parameters such as the Sersic index
of the bulge, the disk scale length, and the disk, bulge, and halo masses. We
also calculate the probability distribution function of the local dark matter
velocity dispersion and density, two quantities of paramount significance for
terrestrial dark matter detection experiments.
  Though the Milky Way models in our chain all satisfy the prescribed
observational constraints, they vary considerably in key structural parameters
and therefore respond differently to non-axisymmetric perturbations. We
simulate the evolution of twenty-five models which have different Toomre Q and
Goldreich-Tremaine X parameters. Virtually all of these models form a bar,
though some, more quickly than others. The bar pattern speeds are ~ 40 - 50
km/s/kpc at the time when they form and then decrease, presumably due to
coupling of the bar with the halo. Since the Galactic bar has a pattern speed
~50 km/s/kpc we conclude that it must have formed recently.
