  The data obtained by the recent modern sky surveys enable detailed studies of
the stellar distribution in the multi-dimensional space spanned by spatial
coordinates, velocity and metallicity, from the solar neighborhood all the way
out to the outer Milky Way halo. While these results represent exciting
observational breakthroughs, their interpretation is not simple. For example,
traditional decomposition of the thin and thick disks predicts a strong
correlation in metallicity and kinematics at $\sim$1 kpc from the Galactic
plane; however, recent SDSS--based work has demonstrated an absence of this
correlation for disk stars. Instead, the variation of the metallicity and
rotational velocity distributions can be modeled using non--Gaussian functions
that retain their shapes and only shift as the distance from the mid--plane
increases. To fully contextualize these recent observational results, a
detailed comparison with sophisticated numerical models is necessary. Modern
simulations have sufficient resolution and physical detail to study the
formation of stellar disks and spheroids over a large baseline of masses and
cosmic ages. We discuss preliminary comparisons of various observed maps and
N--body model predictions and find them encouraging. In particular, the N--body
disk models of Ro\v{s}kar et al. \cite{Roskar 2008} reproduce a change of disk
scale height reminiscent of thin/thick disk decomposition, as well as
metallicity and rotational velocity gradients, while not inducing a correlation
of the latter two quantities, in qualitative agreement with SDSS observations.
