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 super(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. Although 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 nonaxisymmetric perturbations. We simulate the evolution of 25 models which have different Toomre Q and Goldreich- Tremaine X parameters. Virtually all of these models form a bar, although some more quickly than others. The bar pattern speeds are image40-50 km s super(-1) kpc super(-1) 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 image50 km s super(-1) kpc super(-1), we conclude that it must have formed recently.