  Pulsed emission from almost one hundred millisecond pulsars (MSPs) has been
detected in $\gamma$-rays by the Fermi Large-Area Telescope. The global
properties of this population remain relatively unconstrained despite many
attempts to model their spatial and luminosity distributions. We perform here a
self-consistent Bayesian analysis of both the spatial distribution and
luminosity function simultaneously. Distance uncertainties, arising from errors
in the parallax measurement or Galactic electron-density model, are
marginalized over. We provide a public Python package for calculating distance
uncertainties to pulsars derived using the dispersion measure by accounting for
the uncertainties in Galactic electron-density model YMW16. Finally, we use
multiple parameterizations for the MSP population and perform Bayesian model
comparison, finding that a broken power law luminosity function with Lorimer
spatial profile are preferred over multiple other parameterizations used in the
past. The best-fit spatial distribution and number of $\gamma$-ray MSPs is
consistent with results for the radio population of MSPs.
