\section{Introduction}
In a universe devoid of measurement errors, estimating distances from parallaxes would involve simply taking the reciprocal of the measured parallax. However, we do not live in such universe. Measurement errors are always present and estimates must take these into account. This is the main question of this paper: Given a measured parallax $\varpi$ and its measurement uncertainty $\sigma_\varpi$, how can we best estimate the distance $r$ and its uncertainty?
\begin{figure*}
\includegraphics[width=\hsize]{plot_varpi_and_r_distributions_all.pdf}
\caption{The left panel is the probability $P(\varpi|r,\sigma_\varpi)$ as a function of $\varpi$. The middle panel is $P(\varpi|r,\sigma_\varpi)$ as a function of $r$ for an object with $f_{\rm obs}=0.2$, and the right panel is as the middle panel but for $f_{\rm obs}=0.5$. The shaded areas indicate the $2\sigma$ credible interval about $1/r$ and the corresponding transformed credible interval.}
\label{fig:rDist}
\end{figure*}
\begin{figure}
\includegraphics[width=\hsize]{plot_GUMS_fObsDistributionLinear.pdf}
\caption{The probability density (top row) and cumulative distribution (bottom row) of $f_{\rm obs}$ (solid lines) and $f\true$ (dashed lines) of stars in the \textit{Hipparcos}{} (blue lines) and the GUMS (red lines) catalogue. $f\true$ is strictly positive by definition. Note that the plot only covers a subrange of all possible $f$ values. GUMS \citep{rob12} is a synthetic catalogue that simulates the expected content of the \textit{Gaia}{} catalogue, and will be discussed further in Sect.~\ref{subsec:gums}.}
\label{fig:fObsDist}
\end{figure}

This question is of pressing importance, because upcoming data releases from the \textit{Gaia}{} mission will measure parallaxes for all objects down to its magnitude limit of about $G=20$ (approximately $10^9$ objects) with expected accuracies of 0.01\,mas--1\,mas \citep{lin12}. The Large Synoptic Survey Telescope (LSST) hopes to extend this by 4 magnitude with a parallax accuracy of the order of 1\,mas \citep{ive08}. 

A common practice is to report $1/\varpi$ as the distance estimate and give its uncertainty from a first-order Taylor expansion as $\sigma_\varpi/\varpi^2$. But this is problematic because $\varpi$ is noisy and the transformation to distance is highly nonlinear. If we assume a Gaussian noise model for our parallaxes, then the measurement can be considered to have been drawn from a Gaussian with mean $1/r$ and standard deviation $\sigma_\varpi$
\begin{equation}
P(\varpi|r,\sigma_\varpi) = \frac{1}{\sqrt{2\pi}\sigma_\varpi}\exp\left[-\frac{1}{2\sigma_\varpi^2}\left(\varpi-\frac{1}{r}\right)^2\right],\;\sigma_\varpi\geq 0.
\label{eq:likelihood}
\end{equation}
We see in the left panel of Fig.~\ref{fig:rDist} that while this distribution is symmetric in $\varpi$, it is not in $r$. The distribution is skewed even for an object with a small fractional parallax error $f_{\rm obs}=\sigma_\varpi/\varpi$ of 0.2 (central panel of Fig.~\ref{fig:rDist}). Because the likelihood does not go to zero even at infinite distance, confidence intervals could extend up to infinity even for finite parallax errors, giving meaningless results (see \cite{cbj15} for more discussion). This is illustrated in the right-panel of Fig.~\ref{fig:rDist}. For such skewed distributions it is not obvious that the mode is even a good summary.

Another problem with using $1/\varpi$ as the distance estimator is that it does not work for parallaxes with negative values. 
The measurement model of Eq.~\ref{eq:likelihood} has a nonzero probability of drawing negative parallaxes, and this probability grows with increasing fractional parallax error and increasing distance.
Such measurements occur in reality, and tell us that the source is distant and/or has a large measurement error. 

One way to ameliorate these problems is by discarding negative parallaxes and retaining only stars with high accuracy (say $f_{\rm obs}\leq 0.2$). However, as we can see from Fig.~\ref{fig:fObsDist} (bottom panel, solid red line), we expect 80\% of objects in the \textit{Gaia}{} catalogue to have either negative parallaxes or $f_{\rm obs}>0.2$. This is 
quite different from the \textit{Hipparcos}{} catalogue \citep{per97}, in which about 55\% of the sources have $f_{\rm obs}\geq 0.2$ (bottom panel of Fig.~\ref{fig:fObsDist}, solid blue line). Discarding data with large measurement errors is a waste of hard-won data, and will bias subsequent analyses by removing fainter or farther sources \citep{are99, smi06}.

There are of course instances when one is not interested in distances explicitly. Sometimes we want to fit a model and can often do this by using it to predict parallaxes which we then compare directly with the measured parallaxes (e.g. \cite{rom52, jun71, van98}). In this case working in the parallax-domain is preferable as it is much easier to deal with the uncertainties. But here we have no model and are interested specifically in the problem of distance inference for single objects.

In a previous tutorial paper, one of us looked at this problem and showed that one has to take a proper inference approach which cannot avoid the specification of a prior \citep[henceforth Paper I]{cbj15}. This paper examined three priors on $r$: the uniform distance{} prior, the uniform space density{} prior, and the exponentially decreasing space density{} prior. The properties and performance of distance estimators based on these priors was then assessed using noisy toy data. This showed the problems of using naive priors like the uniform distance{} prior and suggested instead that priors which converge asymptotically to zero at infinite distance are to be preferred.

Here we extend this analysis in two important ways. First, we do distance inference with a number of priors and estimators on a simulation of the entire \textit{Gaia}{} catalogue (Sect.~\ref{subsec:gums}--\ref{subsec:simu}). Second, we introduce a new prior based on our knowledge of the distribution of stars in the Galaxy which includes their observability (including accounting for extinction effects) by \textit{Gaia}{} (Sect.~\ref{subsec:mw}). We assess the performance of the four priors in a manner similar to Paper I (Sect.~\ref{sec:results}), and further investigate the reliability of distance uncertainties based on credible intervals derived from the posterior distributions (Sect.~\ref{subsec:formalErrors}). We also considers the inclusion of photometric measurements and a Hertzsprung-Russell Diagram (HRD) prior together with the Milky Way{} prior and compare its performance to parallax-only distance inference (Sect.~\ref{sec:mwhrd}).

A summary of the notations used in this paper is shown in Table~\ref{tab:notation}.
\begin{table}
\caption{Notations used in this paper.}
\label{tab:notation}
\begin{center}
\begin{tabular}{ll}
\hline\hline
$\varpi$ & Measured parallax\\
$\sigma_\varpi$ & Parallax measurement uncertainty\\
$r\true$ & True distance\\
$r\median$ & Distance from the median of the posterior\\
$r\mode$ & Distance from the mode of the posterior\\
$r_{\rm lim}$ & Limiting distance imposed in several priors\\
$f_{\rm obs}$ & $\sigma_\varpi/\varpi$, observed fractional parallax error\\
$f\true$ & $\errVarpir\true$, true fractional parallax error\\
$P(\theta)$ & Probability distribution function (PDF) of $\theta$\\
$P(\theta|x)$ & Posterior distribution of $\theta$ given $x$\\
$P^*(\theta)$ & Unnormalised PDF\\
$C(\theta)$ & Cumulative distribution function (CDF) of $\theta$\\
\hline
\end{tabular}
\end{center}
\end{table}

\section{The methods and the models}
\label{sec:method}
\subsection{Distance inference}
Suppose that a star at true distance $r$ is measured to have a parallax $\varpi$ with measurement error $\sigma_\varpi$. We want to find the posterior probability density (PDF) $P(r|\varpi,\sigma_\varpi)$. Using Bayes' theorem this is expressed in terms of the likelihood of the observation $P(\varpi,\sigma_\varpi)$, and the prior probability $P(r)$ of the distance, as
\begin{equation}
P(r|\varpi,\sigma_\varpi) = \frac{1}{Z}P(\varpi|r,\sigma_\varpi)P(r),
\end{equation}
where
\begin{equation}
Z = \int^\infty_0 dr P(\varpi|r,\sigma_\varpi)P(r)
\end{equation}
is the normalisation constant.

The likelihood is given in Eq.~\ref{eq:likelihood}. This expresses the probability density over the measurement---the parallax---given the true distance, and is the one used in the \textit{Gaia}{} data processing.

The prior expresses our knowledge and assumptions about the distance distribution, independent of the specific measurement. As was explained in Paper I, we have essential knowledge which should be used if parallaxes are to be interpreted correctly in the general case.

\begin{table*}
\caption{The three isotropic priors used in this paper, summarised from \cite{cbj15}.}
\label{tab:priors}
\begin{center}
\begin{tabular}{>{\raggedright\arraybackslash}p{3.6cm}ll}
\hline\hline
Prior & Form & Mode of the posterior\\
\hline
Uniform distance & 
$P_{\rm ud}(r) = 
\begin{dcases*}
\frac{1}{r_{\rm lim}} & for $0<r\leq r_{\rm lim}$,\\
0 & otherwise.
\end{dcases*}$
&
$r\mode = 
\begin{dcases*}
\frac{1}{\varpi} & for $0<\frac{1}{\varpi}\leqr_{\rm lim}$,\\
r_{\rm lim} & for $\frac{1}{\varpi}>r_{\rm lim}$ or $\varpi\leq 0$.\\
\end{dcases*}$
\\
\hline
Uniform space density & 
$P_{\rm usd}(r) = 
\begin{dcases*}
\frac{3}{r_{\rm lim}^3}r^{2} & for $0<r\leq r_{\rm lim}$,\\
0 & otherwise.
\end{dcases*}$ &
$\begin{aligned}
r_{\rm mode} &= \frac{1}{4f_{\rm obs}^2\varpi}\left(1-\sqrt{1-8f_{\rm obs}^2}\right)\\
r\mode &=
\begin{dcases*}
r_{\rm mode} & for $\varpi>0$, $f_{\rm obs}<1/\sqrt{8}$, $r_{\rm mode}\leqr_{\rm lim}$,\\
r_{\rm lim} & otherwise.
\end{dcases*}
\end{aligned}$
\\
\hline
Exponentially decreasing space density &
$P_{\rm exp}(r) =
\begin{dcases*}
\frac{1}{2L^3}r^2\exp\left(-\frac{r}{L}\right) & for $r>0$,\\
0 & otherwise.
\end{dcases*}$
&
Solve $\frac{r^3}{L}-2r^2+\frac{\varpi}{\sigma_\varpi^2}r-\frac{1}{\sigma_\varpi^2} = 0$ for $r$.
\\
\hline
\end{tabular}
\end{center}
\end{table*}
\begin{figure*}
\includegraphics[width=\hsize]{plot_posteriors_wCDF.pdf}
\caption{The unnormalised posterior PDF (top row) and the corresponding CDF (bottom row) for the uniform distance{} prior (left column), uniform space density{} prior (middle column), and the exponentially decreasing space density{} prior (right column). The posteriors are shown for various $f_{\rm obs}$ as indicated by the value near the corresponding line. Positive $f_{\rm obs}$ are shown in black lines, while negative $f_{\rm obs}$ are shown in red lines. Green dashed lines indicate the corresponding priors, but $P^*_{\rm ud}(r)$ is not shown because it is simply a uniform function in $r$. The magnitude of the observed parallax is $|\varpi| = 10$\,mas in all cases. For the uniform distance{} and the uniform space density{} priors, the distance cutoff is $r_{\rm lim}=1$\,kpc and for the exponentially decreasing space density{} prior the scale length is $L=1$\,kpc. Note that the PDF and CDF of the exponentially decreasing space density{} prior are shown on a logarithmic scale in $r$.}
\label{fig:posteriors}
\end{figure*}
\subsection{The three isotropic priors}
\label{subsec:3priors}
The first three priors considered in this paper are the three isotropic priors discussed in Paper I. These are the uniform distance{} prior, the uniform space density{} prior, and the exponentially decreasing space density{} prior. We will introduce the anisotropic Milky Way{} prior in Sect.~\ref{subsec:mw}.

The form of the priors and expression for their mode is given is Table~\ref{tab:priors}, and plots for various values of $f_{\rm obs}$ are shown in Fig.~\ref{fig:posteriors}. Since these priors have been discussed extensively in Paper I, readers are referred to that paper for details. We just note here some particular behaviours.

Using the uniform distance{} prior, the mode of the posterior is always $r\mode = 1/\varpi$, except for negative parallaxes or when $1/\varpi>r_{\rm lim}$, in which case $r\mode=r_{\rm lim}$. This is markedly different from the other priors, in which $r\mode$ is also influenced by the parallax uncertainty. We can see this in the top rows of Fig.~\ref{fig:posteriors}. The mode of the posteriors with the uniform space density{} and the exponentially decreasing space density{} priors shift to larger distances as $f_{\rm obs}$ increases, while the mode of the posteriors with the uniform distance{} priors are independent of $f_{\rm obs}$.

The median distance of all the posteriors, on the other hand, shifts to larger distances as $f_{\rm obs}$ increases. The median distance $r\median$ of the posterior with the uniform space density{} prior shifts faster than the posterior using the uniform distance{} because the uniform space density{} prior assumes that there are more objects located at large distances.

In place of the sharp cut-off at $r=r_{\rm lim}$, the exponentially decreasing space density{} prior drops asymptotically to zero as $r\rightarrow\infty$. While more desirable, this can produce a bimodality over a small range of $f_{\rm obs}$ as a transition occurs from ``data-dominated'' posterior to a ``prior-dominated'' posterior. In the right column of Fig.~\ref{fig:posteriors}, this occurs at $f_{\rm obs}\sim 0.3$. When $f_{\rm obs}$ gets large enough the data are less informative, so the bimodality disappears and the posterior is dominated by the prior, as can be seen for the case of $f_{\rm obs}=0.5$ and $f_{\rm obs}=1$. For $f\gtrsim 1$ the posterior becomes indistinguishable from the prior (shown in the green line). This is just the behaviour we expect and want.

A negative parallax is an indicator that the object is located at a large but highly uncertain distance. Thus from the assumed prior we can still obtain information on the distance, and our confidence in this distance estimate depends on $f_{\rm obs}$. In Fig.~\ref{fig:posteriors} (red lines) we show a negative parallax of $\varpi=-10$\,mas with $f_{\rm obs}=-0.1$ and $f_{\rm obs}=-0.25$.
As the negative $f_{\rm obs}$ shifts toward zero, the mode of the posterior shift toward larger distances. This is because as $|f_{\rm obs}|$ gets smaller, we are more confident that the true parallax is close to zero and the object is very far away. As $|f_{\rm obs}|$ gets larger, the mode shifts to smaller distances and converges on the prior.

\begin{figure*}
\includegraphics[width=\hsize]{plot_GUMS_pObs_ingredient.pdf}
\caption{Left: The smoothed intrinsic luminosity function $\phi(M_G)$ (black line, left scale) used in the Milky Way{} prior, and its cumulative distribution function $\Phi(M_G)$ (blue line, right scale). The luminosity function is smoothed using kernel density estimation. Middle: The average line-of-sight $V$-band extinction as a function of distance $r$, shown as the solid black line. The shaded area indicates the full range of extinctions and the dashed lines show the 5\% and 95\% quantiles of the distribution. Right: The line-of-sight fraction of observable stars $p_{\rm obs}(r)$ as a function of distance. The solid black line indicates $p_{\rm obs}(r)$ when extinction effects are neglected, while the shaded area indicates all the possible shape of $p_{\rm obs}(r)$ for all possible line-of-sights. The dashed line show the lower 5\% of $p_{\rm obs}(r)$ (i.e.\ 95\% of all stars at a given $r$ are above the line).}
\label{fig:pObs}
\end{figure*}
\begin{figure*}
\includegraphics[width=\hsize]{plot_MWSpaceDensityWithExtMapPrior_linearScaleY_new.pdf}
\caption{The Milky Way{} prior $P_{\rm MW}^*(r,l,b)$ along the LOS towards the $(l,b)$ directions written on the top and side of the grid, here as a function of heliocentric distance $r$. Each individual components are drawn in different colors, while the black line indicates the total density. The prior is scaled so that the peak distribution is unity.}
\label{fig:prior_mw}
\end{figure*}
\begin{figure*}
\includegraphics[width=\hsize]{{plot_posterior_04_MilkyWayDensityModelPriorWithExtinctionMap_varpi1.000000_newPDF}.pdf}
\caption{The unnormalised posterior PDF $P_{\rm MW}^*(r|\varpi,\sigma_\varpi,l,b)$ for various $f_{\rm obs}$ as indicated by the legend, and various LOS along the $(l,b)$ direction written on the top and side of the grid. Posteriors with positive $f_{\rm obs}$ are drawn in various shades of blue, while negative $f_{\rm obs}$ are drawn in red. Here the observed parallaxes are all kept the same at $|\varpi|=1$\,mas. The black line shows the prior $P_{\rm MW}(r,l,b)$ previously shown in Fig.~\ref{fig:prior_mw}. The PDF are scaled such that they peak at unity. We can see here that the transition of the posterior to the prior for large $f_{\rm obs}$ also happens, but the shape of the posterior depends also on the LOS direction. Note that the narrow peak at $(l,b) = (0\degr,0\degr)$ contains multiple lines.}
\label{fig:post_mw}
\end{figure*}
\begin{figure*}
\includegraphics[width=\hsize]{{plot_posterior_04_MilkyWayDensityModelPriorWithExtinctionMap_varpi0.100000_newPDF}.pdf}
\caption{As Fig.~\ref{fig:post_mw}, but for $|\varpi|=0.1$\,mas. Note that the range of the horizontal axes is different than those in Fig.~\ref{fig:post_mw}.}
\label{fig:post_mw2}
\end{figure*}

\subsection{The Milky Way{} prior}\label{subsec:mw}
The above priors have the advantage of being simple, but the disadvantage of being very different from what we already know about the distribution of stars in the Milky Way.
We introduce here a prior based on a three-dimensional density model $\rho_{\rm MW}(\mathbf{r})$ of the Milky Way which also takes into account selection effects of a magnitude-limited survey such as \textit{Gaia}{}. To model these selection effects we adopt a universal luminosity function $\phi(M_G)$ and an extinction map. 
Since now we include directional information in addition to the measured parallax, we can hope that this will improve our distance estimation.

In constructing the Milky Way{} prior, we assume that the Milky Way has a stellar number density (the number of stars in a volume element $dV$) described by $\rho_{\rm MW}(r,l,b)$, and a universal luminosity function (the fraction of stars with absolute magnitude between $M_G$ and $M_G+dM_G$) described by $\phi(M_G)$. The probability density of observable stars in a volume element $dV$ is then
\begin{equation}
P^*_{\rm MW}(V) = \rho_{\rm MW}(r,l,b) p_{\rm obs}(r,l,b),
\end{equation}
where
\begin{equation}
p_{\rm obs}(r,l,b) = \int^{M_{G,{\rm faint}}(r,l,b)}_{M_{G,{\rm bright}}} dM_G\phi(M_G)
\label{eq:pObs}
\end{equation}
describes the fraction of observable stars along the line of sight (LOS) direction $(l,b)$ as a function of distance $r$.
Here $M_{G,{\rm bright}}$ is the bright end of the luminosity function, and
\begin{equation}
M_{G,{\rm faint}}(r,l,b) = m_{G,{\rm lim}} - 5\log_{10}r + 5 - A_G(r,l,b)
\label{eq:faintlimit}
\end{equation}
is the faint observable limit of the luminosity function, which is a function of distance $r$ and direction $(l,b)$, and depends on both the limiting magnitude $m_{G,{\rm lim}}$ of \textit{Gaia}{} (which we set to 20) and the $G$-band extinction $A_G$. 
We could neglect extinction effects by setting $A_G = 0$, thus making $M_{\rm lim}$ the faintest possible limit. But as we shall see later, this overestimates the limit at low galactic latitudes, especially in regions close to the galactic centre. Neglecting extinction effects here would erroneously increase our visibility towards the galactic centre and consequently risk providing an unrealistic distance inference for a large number of stars. We therefore include extinction effects in our prior. The full form of the Milky Way{} prior is then
\begin{equation}
P^*_{\rm MW}(r, l, b) = r^2\rho_{\rm MW}(r,l,b)p_{\rm obs}(r,l,b).
\label{eq:prior_mw}
\end{equation}
The stellar number density $\rho_{\rm MW}(r,l,b)$ can then be expressed further as the sum of three components: the bulge $\rho_b$, the disc $\rho_d$, and the halo $\rho_h$. These are detailed in Appendix~\ref{app:mwmodel}.

Earlier we defined $\rho_{\rm MW}(r,l,b)$ as the stellar number density of the Milky Way. However, the model that we use actually describes the stellar mass density. We are therefore assuming that the stellar mass traces the stellar counts in the same way everywhere. We are not assuming that all stars have the same mass.

To construct the fraction of observable stars $p_{\rm obs}(r,l,b)$, we use the the stellar colour-magnitude diagram (CMD)\footnote{We only use the absolute magnitude distribution in the GUMS catalogue for the construction of our Milky Way{} prior. We do not use the spatial distribution} from \cite{rob12}, which is modelled using the initial mass function (IMF) and star formation rate (SFR) observed in the solar neighborhood. The CMD for the thin disc assumes multiple formation epochs with its own IMF and SFR history, while the bulge, thick disc, and halo population assumes a single burst of star formation (the GUMS catalogue and its construction will be described in more detail in Sect.~\ref{subsec:gums}). In this paper, however, we combine all these diagrams into a single luminosity function $\phi(M_G)$.

We transform the $V$-band absolute magnitude $M_V$ into the $G$-band absolute magnitude $M_G$ using the transformation polynomial from \cite{jor10}:
\begin{multline}
M_G = M_V - 0.0257 - 0.0924(V-I_C)\\- 0.1623(V-I_C)^2 + 0.0090(V-I_C)^3 \ .
\end{multline}
To construct a smooth luminosity function from the values extracted from the binned CMD, 
we use kernel density estimation weighted by the number of stars per bin.
We use a Gaussian kernel with kernel width $\sigma = 0.25$ magnitudes. The resulting luminosity function $\phi(M_G)$ and its CDF $\Phi(M_G)$ is shown in the left panel of Fig.~\ref{fig:pObs}. The black line shows the smoothed function, while the blue line shows the CDF of the luminosity function. We can then use this CDF with Eq.~\ref{eq:pObs} to calculate the fraction of observable stars $p_{\rm obs}(r,l,b)$, but in order to calculate the faint end of the integration, in Eq.~\ref{eq:faintlimit} we have to take extinction effects into account. For this we use the extinction map of \cite{dri03}, examples of which are shown in Fig.~\ref{fig:extMap}.

The central panel of Fig.~\ref{fig:pObs} shows the LOS extinction as a function of distance $r$, averaged over all directions. Since the variation in $A_V$ can cover approximately 4 orders of magnitude, it is thus important to take extinction into account. Using the extinction map we can determine $M_{G,{\rm faint}}(r,l,b)$ for any position in the Galaxy. The $V$-band extinction $A_V$ is converted to the $G$-band extinction $A_G$ using the ratio $A_G/A_V=0.695$, a factor we computed from simulations of stellar spectra redenned using the \cite{fit99} extinction curve.

We now perform the integral in Eq.~\ref{eq:pObs} to get $p_{\rm obs}(r,l,b)$. This is shown in 
the right panel of Fig.~\ref{fig:pObs} as a function of $r$. The shaded area shows the full range of observabilities across all $(l,b)$. At 1\,kpc, for example, we can see between about 15\% and 50\% of all stars in the GUMS simulation, depending on direction.

Combining $\rho_{\rm MW}(r,l,b)$ and $p_{\rm obs}(r,l,b)$ we can then calculate the Milky Way{} prior $P^*_{\rm MW}(r,l,b)$ using Eq.~\ref{eq:prior_mw}. This is shown in Fig.~\ref{fig:prior_mw} for 16 LOS directions in galactic coordinates $(l,b)$. Due to symmetry in the model of the mass distribution, we only show the LOS directions in the northeastern quadrant of the Milky Way. The variations in the other quadrants will be similar, barring some difference in details due to the extinction. The prior shows a direction dependence, but it is clear that the (thin) disc is the most dominant component in almost any direction, being the most massive among the others. At higher latitudes and distances, the halo starts to be discernible, as does the thick disc at intermediate distances. The bottom right panel shows the prior towards the galactic centre. We see a very large and narrow peak around 8\,kpc, which is the galactic centre. The peak is mostly dominated by the bulge. Had we not included extinction in our prior, this peak would be even larger due to the enormous density of stars in that direction, which would led to erroneous features in the posteriors. 

Example posteriors resulting from these priors are shown in Fig.~\ref{fig:post_mw}, for the same LOS directions as Fig.~\ref{fig:prior_mw}. The measured parallax is $|\varpi|=1$\,mas in all cases. The different lines in each panel show different values of $f_{\rm obs}$.
Here the posterior is almost always unimodal, but nevertheless bimodalities can still appear, for example at $(l,b) = (0\degr, 30\degr)$. The bimodalities are a consequence of both the prior and the likelihood (we saw bimodalities with the exponentially decreasing space density{} which was unimodal).
In general bimodalities only appear when $f_{\rm obs}\gtrsim 0.5$, and only in the directions where there are comparable contributions from multiple components.
If we look at the previously mentioned LOS direction $(l,b) = (0\degr,30\degr)$ again, at $f_{\rm obs}=1$ a second peak in the posterior appears, which indicates that---given this very large measurement error---the posterior considers the possibility that the source is a thick disc object.

Fig.~\ref{fig:post_mw2} shows the same plots but with a different measured parallax, $\varpi = 0.1$\,mas. Here the posteriors are always unimodal.

\begin{figure}
\includegraphics[width=\hsize]{plot_GUMS_spatialDistribution_healpix_nSide0512_logScale_alternate.pdf}
\caption{The Mollweide projection of the spatial distribution of all stellar objects in the GUMS catalogue with $G\leq 20$, in the galactic coordinate system.}
\label{fig:gums}
\end{figure}

\begin{figure}
\includegraphics[width=\hsize]{plot_GUMS_trueDistanceDistribution_allDirections_1plot.pdf}
\caption{The distribution of the true distances $r\true$ of all stellar objects in the GUMS catalogue (black line). The other colored lines indicates the exponentially decreasing space density{} prior for different scale lengths $L$. The PDFs are scaled such that they peak at unity.}
\label{fig:gums_rdist}
\end{figure}
\begin{figure*}
\includegraphics[width=\hsize]{plot_GUMS_trueDistanceDistribution_selectedDirections_linearScaleY.pdf}
\caption{The distribution of the true distances $r\true$ of all stellar objects in the GUMS catalogue, for selected LOS directions (as used in Figs.~\ref{fig:prior_mw}--\ref{fig:post_mw2}. Each LOS for $b>0$ comprises stars in an area of 13.43 square degrees (all stars at \texttt{healpix} \citep{gor05} level $n_{\rm side} = 16$ centered on that direction). At $b=0$ we use $n_{\rm side} = 128$ which corresponds to an area of 0.21 square degrees. The black lines show the Milky Way{} prior for comparison. Both distributions are scaled such that they peak at unity in each panel. The sudden drop in the GUMS distribution at $b=0^\circ$ and $l=90^\circ,\ 135^\circ$, and 180$^\circ$ is an artefact of the GUMS catalogue, which terminates the disc at a galactocentric distance of 14\,kpc \citep{rob12}.}
\label{fig:gums_rdist_directional}
\end{figure*}

\subsection{The simulated data: The GUMS catalogue}
\label{subsec:gums}
We want to investigate the performance of the priors described in Sect.~\ref{subsec:3priors}--\ref{subsec:mw} on objects in a \textit{Gaia}{}-like catalogue. For this we can use the \textit{Gaia}{} Universe Model Snapshot (GUMS) catalogue \citep{rob12}, which simulates what \textit{Gaia}{} was expected to observe (prior to launch) down to its limiting magnitude.

Stars are generated in GUMS using the Besan\c{c}on galaxy model (BGM, \citealt{rob03}). It is based on a four component model of the Galaxy using an evolutionary model, with a specified initial mass function and star formation history. Interstellar extinction is calculated from the map of \cite{dri03}. For each star GUMS reports the apparent magnitudes, colours, positions and distance, as well as other quantities not required in our present work. GUMS only simulates objects that can be reasonably well-observed by \textit{Gaia}{}, so it omits objects fainter than $G=20$ or more than 100\,kpc from the Sun. The final catalogue contains $\sim$$1.6\times 10^9$ stellar objects, the spatial distribution of which are shown in Fig.~\ref{fig:gums}.

Fig.~\ref{fig:gums_rdist} show the distribution of the true distances $r\true$ of all stellar objects in the GUMS catalogue. The distribution peaks at $\sim$2.7\,kpc, which corresponds to the scale length of the thin disc adopted by \cite{rob12}. A small peak is also present at $\sim$6\,kpc, which corresponds to the maximum distribution of observable stars towards the galactic centre (the galactic centre itself is invisible due to extinction).
We overplot the exponentially decreasing space density{} prior for 4 different scale length $L$, $L=\{0.1, 0.5, 1.35, 5\}$\,kpc. None of these distributions closely resembles the GUMS distribution, but the mode of the prior for $L=1.35$\,kpc peaks at the same point as those of the GUMS distribution. 

Fig.~\ref{fig:gums_rdist_directional} shows the distance distribution in the GUMS model for different directions (the same as previously used in Figs.~\ref{fig:prior_mw}--\ref{fig:post_mw2}), which we also compare to our Milky Way{} prior. Being a prior, of course we do not expect it to agree exactly with (our simulation of) reality! Yet it is realistic to believe that our current knowledge can produces a distance distribution which is in broad agreement with what \textit{Gaia}{} will observe. Overall we see a good match between GUMS and our model, but nonetheless some deviations can be seen. The most significant deviations occur at low latitudes in directions toward the galactic plane at intermediate and large distances, but in other directions we see nevertheless small devations at large distances. This is due to the different adopted model for the thick disc and the halo (and in the case of $(l,b)=(0\degr,0\degr)$, the difference is due to the different bulge-to-disc mass ratio), and the fact that we are using a single universal HRD instead of an HRD for each component of the Milky Way. Such prior mismatches are bound to occur in reality, so using GUMS will be a meaningful test of the performance of this prior. 

\begin{figure}
\includegraphics[width=\hsize]{plot_nTransit_5yr_galactic_paper.pdf}
\caption{The number of transits as a function of $(l,b)$ assuming a nominal \textit{Gaia}{} scanning law for five years. Data kindly provided by Berry Holl.}
\label{fig:ntr}
\end{figure}

\begin{figure}
\includegraphics[width=\hsize]{plot_GUMS_08b_sigmaVarpi_GMagRTrue_2dHist.pdf}
\caption{The distribution of the parallax measurement errors $\sigma_\varpi$ for simulated \textit{Gaia}{} observations of the GUMS catalogue. 
The colour shows the probability density per unit log parallax error and per unit magnitude (left) or per unit distance (right).
The diagonal lines in the right panel show the locii of constant $f\true = \errVarpir\true$ as indicated by the labels.}
\label{fig:errVarpi_GR}
\end{figure}

\begin{figure}
\includegraphics[width=\hsize]{plot_GUMS_07c_varpi_sigmaVarpi_2dHist_logScale_logBin_positiveNegative.pdf}
\caption{The distribution of $\sigma_\varpi$ as shown in Fig.~\ref{fig:errVarpi_GR},
but now as a function of the simulated measured parallaxes, $\varpi$.
The colour shows the probability density per unit log parallax error and per unit log parallax.
The left panel is for negative parallaxes while the right panel is for positive parallaxes. 
The diagonal lines show the locii of constant $f_{\rm obs} = \sigma_\varpi/\varpi$ as indicated by the labels.
For negative parallaxes, no stars have measurement errors less than 0.01\,mas.}
\label{fig:varpiSigmaVarpiDist}
\end{figure}

\subsection{The parallax noise model}
GUMS provides noise-free measurements. For our simulation we need to adopt a noise model for \textit{Gaia}{}. We take into account the increased scattered light levels found after launch, which slightly reduce the astrometric accuracy at the faint end compared to pre-launch expectation.
Using the post-launch noise model from \cite{deb14}, the expected Gaussian one sigma uncertainty
in the parallax ($\sigma_\varpi$ in Eq.~\ref{eq:likelihood})
for a source of magnitude $G$ and colour $(V-I_C)$, observed in 
$n_{\rm tr}$ transits of the \textit{Gaia}{} focal plane is
\begin{multline}
\sigma_\varpi\;[\mu{\rm as}] = \left(\frac{\langle n_{\rm tr}\rangle}{p_{{\rm det},G}(G)n_{\rm tr}}\right)^{1/2}\times\\
 \left(-1.631 + 680.766z + 32.732z^2\right)^{1/2}\times\\ 
 \left[0.986 + (1 - 0.986)(V-I_C)\right].
\label{eq:sigmaVarpi}
\end{multline}
where
\begin{equation}
z = \max\left(10^{0.4(12.09 - 15)}, 10^{0.4(G-15)}\right),
\label{eq:zDeb14}
\end{equation}
and $p_{{\rm det},G}(G)$ is the detection probability per transit (\textit{Gaia}{} does real time detection independently on every transit). For a five year mission, $\langle n_{\rm tr}\rangle$ (the sky-averaged number of transits) is 70, which takes into account mission gaps and dead time. $n_{\rm tr}$ depends on the \textit{Gaia}{} scanning law, and is shown as a function of $(l,b)$ for a five year mission in Fig.~\ref{fig:ntr}. $p_{{\rm det},G}(G)$ drops with fainter magnitudes, but is still 95\% at $G=20$ (Table~10 in \cite{jor10}).

Using the above equation, we calculate $\sigma_\varpi$ for every source in GUMS.
Fig.~\ref{fig:errVarpi_GR} shows the distribution as a function of magnitude and distance.
The spread in $\sigma_\varpi$ for objects with the same $G$ is due primarily to the variable number of transits, 
and to a lesser extent due to colour variations across sources.

At the bright end the curve is largely independent of $G$. 
In reality there is still a dependence due to the activation of Time Delayed Integration (TDI) gates used to avoid saturation. This complication is ignored here and a constant parallax noise floor of $\sigma_\varpi\sim 7$\,$\mu$as is used for sources brighter than $G\leq 12.09$.

In the right panel of Fig.~\ref{fig:errVarpi_GR} we see that almost all stars with $r\true\lesssim 100$\,pc will have $f\true=\errVarpir\true$ below 0.1. The number of stars with $f\true\leq 0.1$ decreases quickly as we move to larger distances.
Only a very small fraction of stars with $r\true\gtrsim 10$\,kpc have $f\true<0.1$. Comparing the
left and right panels, we see that these stars are near the turn-off at $G\sim 12$: they are distant but bright (mostly giants), so relatively accurate parallaxes can still be obtained.

\subsection{Simulated observed catalogue}
\label{subsec:simu}
We now use the data in the GUMS catalogue $(l,b,r\true,G,V-I_C)$ together with the measurement model 
(the likelihood in Eq.~\ref{eq:likelihood}) and the specification of $\sigma_\varpi$ in
Eq.~\ref{eq:sigmaVarpi} and \ref{eq:zDeb14} to simulate the parallax for every star in the GUMS catalogue, by drawing a parallax at random from the likelihood.
Fig.~\ref{fig:varpiSigmaVarpiDist} shows the distribution of the parallax uncertainties again, but now as a function of these measured parallaxes. 
The right panel is for positive parallaxes. There is, of course, a very broad distribution in parallax errors for a given parallax, due to the wide range of actual distances and magnitudes of the stars. Stars with small measured parallaxes tend to be distant and faint, and so generally have larger parallax errors and likewise larger fractional parallax errors. The left panel shows negative parallaxes, which is approximately 15\% of all stars.
We can obtain negative parallaxes when $1/r\true$ is very close to zero and the parallax measurement error $\sigma_\varpi$ is large. We see in the left panel that the minimum $f_{\rm obs}$ for negative parallaxes is around 0.2.

Given the simulated parallaxes and corresponding uncertainties, we compute the posterior PDF over distance for each star, using the four priors described in Sect.~\ref{subsec:3priors} and \ref{subsec:mw}.
For the uniform distance{} and uniform space density{} priors, we set $r_{\rm lim} = 100$\,kpc as the limiting distance. This is a distance large enough that covers all stellar objects within the GUMS catalogue. For the scale length $L$ of the exponentially decreasing space density{}, we take four different values: $L = \{0.1, 0.5, 1.35, 5\}$\,kpc. We already saw in Fig.~\ref{fig:gums_rdist} that for $L=1.35$\,kpc the posterior closely matches the distribution of $r\true$ in the GUMS catalogue, but nevertheless in the case of a mismatch in the scale length we are interested in seeing how it affects the distance determination.

From the posteriors we estimate distances using the mode and median. To calculate the median, as well as a 90\% percent credible interval, the posterior is numerically integrated by sampling on a dense regular grid in $\log r$. The modes for the isotropic priors are given in Table~\ref{tab:priors}, while the mode of the Milky Way{} prior is found using a golden section search on the aforementioned grid, initialised on the point closest to the mode. We do not use the mean, as some posteriors are very skewed, for which it is unrepresentative.

\begin{figure*}
\includegraphics[width=\hsize]{plot_GUMS_comparison_rTrueVSrMo_00_pureAstrometric_logBin.pdf}
\caption{Comparison of the median distance estimate $r\median$ for all stars in the GUMS catalogue to their true distances $r\true$, using the uniform distance{} prior (top left), the uniform space density{} prior (top middle), the Milky Way{} prior (top right), and the exponentially decreasing space density{} prior (four panels at the bottom) with four different values of scale length $L$ as indicated by the legend. The black diagonal line indicates a perfect match between the estimated and true values. The colour shows the probability density per unit log distance squared.}
\label{fig:comparisons_mo}
\end{figure*}
\begin{figure*}
\includegraphics[width=\hsize]{plot_GUMS_comparison_rTrueVSrMo_04_milkyWayPrior_directional_logBin.pdf}
\caption{A comparison between the mode distance $r\mode$ for all stars in the GUMS catalogue with their true distances $r\true$, using the Milky Way{} prior. Here the comparison is drawn for 16 selected LOS directions as previously used. The black lines indicate a perfect match between the estimated and true value. The colour shows the probability density per unit log distance squared. The sudden drop in density at $b=0^\circ$ and $l=90^\circ,\ 135^\circ$, and 180$^\circ$ is due to the termination of the disc in the GUMS catalogue, which is also observed in Fig.~\ref{fig:gums_rdist_directional} (q.v.).}
\label{fig:comparison_mo_mw_dir}
\end{figure*}
\begin{figure}
\includegraphics[width=\hsize]{plot_GUMS_summaryPerformanceFTrue_mode_longRangeF_00_3panels.pdf}
\caption{Performance of the mode distance $r\mode$ for all priors (coloured lines) as a function of $f\true$. The three panels show the variation of the bias $\overline{x}_{Mo}$ (top), the RMS $\overline{x^2}^{1/2}_{Mo}$ (middle), and the standard deviation $\sigma_{x,\mode}$ (bottom).}
\label{fig:summary_performance_ftrue_mode}
\end{figure}
\begin{figure*}
\includegraphics[width=\hsize]{{plot_GUMS_summaryPerformanceFTrue_mode_longRangeF_directional_RMS_00_linearF_magCutoff0.00}.pdf}
\caption{The RMS of the scaled residuals using the mode, $\overline{x^2}^{1/2}_{Mo}$, for the 16 selected LOS directions for all four priors.}
\label{fig:summary_performance_ftrue_mode_dir}
\end{figure*}
\begin{figure}
\includegraphics[width=\hsize]{plot_GUMS_summaryPerformanceFTrue_mode_expDecSpaceDens_3panels.pdf}
\caption{As Fig.~\ref{fig:summary_performance_ftrue_mode}, but for the exponentially decreasing space density{} prior with various scale length $L$. The dark blue line corresponding to $L=1.35$\,kpc is the same one as those in Fig.~\ref{fig:summary_performance_ftrue_mode}. Scale lengths other than $L=1.35$\,kpc generally give poorer performance, as these are a poorer match to the scale lengths in the Galaxy.}
\label{fig:summary_performance_ftrue_mode_exp}
\end{figure}
\section{Test results and discussion}
\label{sec:results}
In this section we analyse the performance of the priors. In the interest of brevity, here we only look at the mode since we found that overall the mode performs better than the median. The results from using the median are presented and discussed in Appendix~\ref{app:comparisons_md}.

\subsection{Comparison with the true values}
Fig.~\ref{fig:comparisons_mo} compares the mode distance estimate $r\mode$ with the true values $r\true$ for all stars, using each of the four prior. As we already saw in Fig.~\ref{fig:errVarpi_GR} (right panel), all stars within 100\,pc will have their distances measured extremely accurately for all priors. Furthermore, there is a high density around the diagonal line (note that the density is on a log scale), indicating that also at large distances a lot of stars achieve accurate distances with all priors. Beyond 10\,kpc, however, this high density around the diagonal line diminishes for all priors, as the number of stars with accurately measured parallaxes get rarer.

The spread in $r\mode$ arises because of the large range of values of $f\true$ and $\varpi$ present at a given value (or rather a small range) of $r\true$. The variation of this spread with $r\true$ is therefore complicated and cannot be easily explained, as it depends on a large number of factors, including the shape of the posterior and the complexity of distribution of stars in the GUMS models.

The distribution of the uniform distance{} prior is of particular interest, because the mode of the posterior PDF with the uniform distance{} prior is the inverse of the observed parallax $\varpi$ (provided it is not too small or negative). Although hard to see in the plot, there are many stars with $r\mode=100$\,kpc, which is the maximum distance $r_{\rm lim}$ imposed by the prior. This is $r\mode$ for negative parallaxes and parallaxes with $1/\varpi>r_{\rm lim}$. We also see the same cluster in the results using the uniform space density{} prior. Here there are more stars with $r\mode=100$\,kpc because not only stars with negative parallaxes have $r\mode=100$\,kpc, but also all stars with $f_{\rm obs}>1/\sqrt{8}$. These aside (which are easily identified in the results), the uniform space density{} prior looks superior to the uniform distance{} prior, because the former has a smaller spread above the $r\mode=r\true$ line.

In the distribution of the Milky Way{} prior we see a horizontal elongation across $r\true$ at around $r\mode\sim 8$\,kpc. These are stars in the direction towards the galactic centre, with very high observed fractional parallax errors $f_{\rm obs}$. For such stars, the posterior is practically the same as the prior (seen in Fig.~\ref{fig:posteriors} and \ref{fig:post_mw}), so the inferred distance is close to the mode of the prior, which in this case is around 8\,kpc (the distance to the galactic centre in this prior). We see some other elongations corresponding to the modes of the prior in other directions, but these are much weaker.

Similar elongation is also observed in the distribution of the exponentially decreasing space density{} prior---shown in the bottom row of Fig.~\ref{fig:comparisons_mo}---at the mode of this prior, $r\mode=2L$.

Since the Milky Way{} prior has different shapes depending on the LOS direction, we show in Fig.~\ref{fig:comparison_mo_mw_dir} the $r\mode$ vs.\ $r\true$ comparison for the same 16 LOS directions used previously (cf. Fig.~\ref{fig:gums_rdist_directional}). Looking at the distribution, it is more clear that there are elongations of stars with poorly measured parallaxes around a preferred mode. At directions away from the galactic center, however, the elongations are closer to the true distance of most stars, indicating that the prior can help in estimating distances of stars with poor parallax measurements, provided that there are reasonably good match between the prior and the underlying true distribution. We have seen that we have a good match between the prior for the thin disc and its corresponding ``true'' distribution in GUMS, and this is advantageous for the prior as a whole since most stars are thin disc objects. In directions around the galactic centre, however, at $(l,b) = (0\degr,0\degr)$ and $(l,b) = (0\degr,30\degr)$, the elongations are relatively more pronounced, which is due to the prior mismatch in this direction. Note the high density around the galactic centre.

\subsection{Statistical performance as a function of $f\true$}
We evaluate the performance of the mode distance estimator $r\mode$ as a function of $f\true$. We divide the data into bins of $f\true$, then for every object, $i$, in the bin we calculate the scaled residual 
\begin{equation}
x_{Mo,i} = \frac{r_{Mo,i} - r_{{\rm true},i}}{r_{{\rm true},i}} \ .
\label{eq:x}
\end{equation}
For each bin we then calculate the bias and standard deviation as explained in Paper I, as well as the root mean square (RMS) of the scaled residual.
The standard deviation measures the scatter about the bias, whereas the RMS measures the entire error including the bias. The variation of these metrics with $f\true$ for the different priors is shown in Fig.~\ref{fig:summary_performance_ftrue_mode}. We can immediately conclude from this plot that both the uniform distance{} and uniform space density{} priors performs badly. The latter, for example, becomes heavily biased even at a low values of $f\true$ ($\sim$0.2), and the errors (in terms of RMS and standard deviation $\sigma_{x,\mode}{}$) rise rapidly too. Contrary to our expectations, however, the Milky Way{} prior does not show the best overall performance; rather it is the exponentially decreasing space density{} prior with $L=1.35$\,kpc. We see that for large $f\true$, the standard deviation of the scaled residual $\sigma_{x,Mo}$ can be as low as $\sim$0.14, and the peak deviation is $\sim$0.46 while for the Milky Way{} prior its peak deviation is $\sim$0.90, and that occurs at a relatively low $f\true$ of $\sim$0.47.

We should be careful in interpreting these results, however. For $f\true\gtrsim 1$, the exponentially decreasing space density{} prior gives a negative bias which increases in magnitude with $f\true$, while the bias in the Milky Way{} prior decreases towards zero. The large negative bias of the former prior arises because of its isotropy and chosen scale length $L$ of 1.35\,kpc, which produces distance estimates for poorly measured parallaxes all at around $2L=2.7$\,kpc. For a majority of stars with large $f\true$, this is far from their true distance. This is the unavoidable consequence of choosing a simple prior. The RMS includes the effect of this bias, and we see in the central panel of Fig.~\ref{fig:summary_performance_ftrue_mode} the RMS of the exponentially decreasing space density{} prior increases with $f\true$ at larger values of $f\true$, whereas the Milky Way{} prior does not.

One must also appreciate that Fig.~\ref{fig:summary_performance_ftrue_mode} is the performance averaged over all directions in the Galaxy, whereas the Milky Way{} prior has a directional dependence. In Fig.~\ref{fig:summary_performance_ftrue_mode_dir} we show the RMS as a function of direction for all priors. For directions away from the galactic centre, the Milky Way{} prior performs better than the exponentially decreasing space density{} for $f\true\lesssim 1$. In contrast, the Milky Way{} prior has worse performance in directions toward the galactic centre. As expected, this is the effect of a mismatch between the bulge model as a prior and the ``true'' bulge used in the simulated observations from GUMS (cf. Fig~\ref{fig:gums_rdist_directional}). We can see that away from the galactic centre, the Milky Way{} prior actually performs best. However, due to the high number of stars in the directions around the galactic centre, the poor performance of the Milky Way{} prior in Fig.~\ref{fig:summary_performance_ftrue_mode} is dominated by stars in these directions. This should be taken into consideration when using the Milky Way{} prior.

In Fig.~\ref{fig:summary_performance_ftrue_mode_dir}, we can also see that the exponentially decreasing space density{} prior with $L=1.35$\,kpc performs best in directions at the galactic plane, as the chosen scale length is a good match with the ``true'' thin disc used in GUMS. Variations in the behaviour of the exponentially decreasing space density{} prior between the panels are due entirely to the different distribution of stars in GUMS, as the prior itself is isotropic.

In terms of standard deviation, the exponentially decreasing space density{} prior performs better than the Milky Way{} prior. In terms of bias the relative performance depends on $f\true$. Overall the exponentially decreasing space density{} prior with $L=1.35$\,kpc is the best performing prior for $f\true\lesssim 2.5$ (which is already a very large error). We show the results of this prior only for $L=1.35$\,kpc because it is the overall best performing (both in averaged and directional performance) among the four scale lengths $L$ investigated, as expected (cf. Fig.~\ref{fig:gums_rdist_directional}. Other scale lengths, however, performs better in different regime of $f\true$. This can be seen in Fig.~\ref{fig:summary_performance_ftrue_mode_exp}. For $f\true\lesssim 0.6$, those with $L=0.5$\,kpc performs best, while for $f\true\gtrsim 2$ it is $L=3$\,kpc.

Note that all of these plots tell us only the performance at a given $f\true$, but not how many stars actually have that value of $f\true$. This can be seen in Fig.~\ref{fig:fObsDist}, and should be taken into mind when comparing performances.

\begin{figure}
\includegraphics[width=\hsize]{plot_GUMS_summaryPerformanceFObs_mode_longRangeF_00_3panels.pdf}
\caption{As Fig.~\ref{fig:summary_performance_ftrue_mode}, but now plotted as a function of $f_{\rm obs}=\sigma_\varpi/\varpi$. Note that now the horizontal axis extends to negative values.}
\label{fig:summary_performance_fobs_mode}
\end{figure}
\begin{figure}
\includegraphics[width=\hsize]{plot_GUMS_06_fTrue_fObs_2dHist_fObsLinear.pdf}
\caption{The distribution of $f_{\rm obs}$ and its corresponding $f\true$. Note that $f\true$ is plotted in log scale, while $f_{\rm obs}$ is in linear scale and extend to negative values. The black line indicates equal values of $f_{\rm obs}$ and $f\true$. The colour shows the probability density per log $f\true$ per $f_{\rm obs}$.}
\label{fig:fObsTrueDist}
\end{figure}
\begin{figure*}
\includegraphics[width=\hsize]{{plot_GUMS_summaryPerformanceFObs_mode_longRangeF_directional_RMS_00_linearF_magCutoff0.00}.pdf}
\caption{As Fig.~\ref{fig:summary_performance_ftrue_mode_dir}, but now plotted against the observed fractional parallax error $f_{\rm obs}=\sigma_\varpi/\varpi$.}
\label{fig:summary_performance_fobs_mode_dir}
\end{figure*}
\begin{figure}
\includegraphics[width=\hsize]{{plot_GUMS_formalErrors_ZHistogram_rMo_00_magCutoff0.00}.pdf}
\caption{The distribution of the standardised residual $z_{_{Mo}}$ (defined in Eq.~\ref{eqn:standardized_residual}), for all priors. The black dashed line indicates the standard normal distribution.} All distributions are scaled to have their maxima at one.
\label{fig:zhist}
\end{figure}
\subsection{Statistical performance as a function of $f_{\rm obs}$}
In the previous section we showed the performance of the distance estimators as a function of $f\true$. These plots are useful for predicting the expected peformance in terms of the {\em expected} fractional parallax errors. But as the true distances are not known in a real application, we will also want to know, for a given {\em measured} fractional parallax error, $f_{\rm obs}=\sigma_\varpi/\varpi$, what the expected performance is (We would still not know the true performance, of course). Note that $\sigma_\varpi$ is not measured as such---we do not know the true uncertainty---but rather is estimated from a noise model.
Fig.~\ref{fig:summary_performance_fobs_mode} shows the performance of the mode now as a function of $f_{\rm obs}$, averaged over all directions. This figure includes negative values of $f_{\rm obs}$ because observed parallaxes can be negative. This is quite a different plot, because the horizontal axis now includes the noisy measurement.

Comparing this to the similar plot against $f\true$ (Fig.~\ref{fig:summary_performance_ftrue_mode}) we instantly see a significant difference: When plotted against $f_{\rm obs}$, the uniform distance{} prior performs much better for positive parallaxes, and in fact as well as or even better than the Milky Way{} prior for $f_{\rm obs}\lesssim 1.2$. The uniform distance{} prior (and uniform space density{} prior) remains useless for negative parallaxes, however, because the mode is then always $r_{\rm lim}=100$\,kpc.

The reason for this change in behaviour for the uniform distance{} prior is two-fold. First, some stars which previously had smaller values of $f\true$ will, due to the noise, achieve larger (absolute) values of $f_{\rm obs}$. Hence the good performance we previously saw at low $f\true$ is partly distributed to larger values of $f_{\rm obs}$
(Noise will also give some stars even lower values of $f_{\rm obs}$, but we already saw good performance at low $f\true$).
This can also happen for the other priors, of course, but these did not previously have such terrible performance at larger values of $f\true$, so the effect is less noticeable. That this does not occur for the uniform space density{} prior is because for $f_{\rm obs}$ larger than about 0.3, the mode is always at $r_{\rm lim}$, which is a very poor estimate for most stars. The uniform distance{} prior, in contrast, at least has a variable mode in those cases ($1/\varpi$), even if it is often a poor and biased estimate (given that $\varpi$ is noisy). 
When we were plotting the performance against $f\true$, the bin can contain stars with negative $f_{\rm obs}$. In Fig.~\ref{fig:fObsTrueDist} we show the distribution of $f_{\rm obs}$ and the corresponding $f\true$. We see that negative $f_{\rm obs}$ can even be drawn already from $f\true\sim 0.2$. This is the second reason for the uniform distance{} prior giving better results when plotted against $f_{\rm obs}$. Stars with negative parallaxes no longer appear in this plot: the line is off to the top left of all the panels in Fig.~\ref{fig:summary_performance_fobs_mode}. 

At negative parallaxes (negative $f_{\rm obs}$), we see that the Milky Way{} prior is less biased than the exponentially decreasing space density{}, but has a larger standard deviation. This is what we might expect, given that in this regime the prior dominates.

As mentioned before, these summary plots do not tell us how many stars there are at each value of $f_{\rm obs}$. This can be seen in Fig.~\ref{fig:fObsDist}.

Fig.~\ref{fig:summary_performance_fobs_mode_dir} shows the directional dependence of the mode estimator now as a function of $f_{\rm obs}$ (cf. Fig.~\ref{fig:summary_performance_ftrue_mode_dir}). The plots are noisy at large $f_{\rm obs}$ for at high latitudes because of the paucity of stars. We now see that the Milky Way{} performs better than any other prior, except for directions toward the galactic centre. This is the same behaviour as in Fig.~\ref{fig:summary_performance_ftrue_mode_dir}.

Overall the exponentially decreasing space density{} prior performs best, but for directions away from the galactic centre the Milky Way{} prior performs best. In general we cannot recommend using the uniform distance{} or uniform space density{} priors.

\subsection{The formal errors}
\label{subsec:formalErrors}
Just as important as achieving accurate distance estimates is estimating accurate uncertainties in those estimates.
An estimate without an error bar is essentially useless. Here we use the 90\% credible interval $[r_{5},r_{95}]$ of the posterior PDF computed symmetrically about the median (i.e.\ the difference between the 95\% and 5\% quantiles). As many people are still used to think about Gaussian-like one sigma errors, we convert this interval into what would be a one sigma standard deviation if the PDF were Gaussian, by scaling the computed 90\% interval to the 68.3\% interval (this also capture the shape of the tailing distributions). We calculate
\begin{equation}
\sigma_r = \frac{\left(r_{95}-r_{5}\right)}{2s}
\label{eqn:standardized_residual}
\end{equation}
where $s = 1.645$ is the ratio of the 90\% to 68.3\% credible interval in a Gaussian. We then calculate the standardised residuals
\begin{equation}
z_{Mo} = \frac{r\mode-r\true}{\sigma_r} \ .
\end{equation}
{\em If} the residuals were Gaussian, we would expect this quantity to have a distribution which is a standardised Gaussian. In any case, for given residuals, values of $z$ much larger than one imply that the uncertainties---as estimated by our credible interval---are underestimated, and values of $z$ much less than one imply that they are overestimated.

Fig.~\ref{fig:zhist} shows the distribution of $z$ for all the priors. The distribution of $z$ for the uniform distance{} and the uniform space density{}
have a component which is centered around zero and very narrow, the latter suggesting that the uncertainties are overestimated. But they are both bimodal with significant peaks at large positive values of $z$. These correspond to the stars with large fractional parallax errors or negative parallaxes, when the mode is pushed to $r_{\rm lim}$. 
We can now see much better how many stars are concentrated into those narrow lines at $r\mode=r_{\rm lim}$
in Fig.~\ref{fig:summary_performance_ftrue_mode_dir} for these two priors.
Thus these priors only achieve unbiased results for a limited fraction of all objects, and even a posteriori we do not know which of the objects with larger $f_{\rm obs}$ these are.

The huge secondary peak of the uniform space density{} prior indicates that the mode is a poor estimator for the uniform space density{} prior.

The distribution for the exponentially decreasing space density{} prior has its mode at $z=0$, but is skewed towards negative $z$.
This is just the negative bias we discussed before. But as the half-width at half-maximum (HWHM) of the distribution is 1.17, which is
close to 1.2 (the value for a Gaussian). At least the uncertainty estimates are reasonable.

The distribution for the Milky Way{} prior is in many ways much closer to what we want for our estimators: centered on zero and symmetric. Its HWHM is 0.84,
which however suggests that the uncertainties might be overestimated. But we must remember that they have been calculated from a 90\% credible interval, so we do not expect the standardised residuals to follow a standardised Gaussian.

Overall, we conclude that although when averaged over all directions the Milky Way{} prior does not perform as well as the exponentially decreasing space density{} prior in terms of distance estimation, it does performs better in terms of uncertainty estimates. This we see from the distribution of the residuals, which centered on zero and is symmetric (see Fig.~\ref{fig:zhist}). The uniform distance{} and uniform space density{} uncertainty estimates are basically useless.

\section{Adding photometric measurements to improve the distance inference}
\label{sec:mwhrd}

Once the fractional parallax error increases above about 100\%, distance estimation with any prior will be inaccurate.
For some types of stars we can get better distance estimates from their colours or spectra, together with photometry, using what are sometimes called photometric distances. This relies on prior information from stellar structure and evolution models, which predict absolute magnitudes based on the colours or spectra. Combined with the measured photometry, and correcting for extinction, a distance estimate can be obtained. For some types of stars, such as red-clumps and giants, distance accuracies of up to respectively 5\%--10\% \citep{bov14} and 20\% \citep{san16} are achieved.

In addition to the $G$-band magnitude, \textit{Gaia}{} will also observe all stars using low resolution spectrophotometry obtained in the blue and red. These will be used to estimate stellar parameters \citep{cbj13} which in turn can give us photometric distances. \textit{Gaia}{} will also provide the integrated fluxes from these to produce the $G_{\rm BP}$ and $G_{\rm RP}$ magnitudes, which are essentally very broad blue and red filters \citep{jor10}. For simplicity we will just examine here how the use of the three magnitudes $G$, $G_{\rm BP}$, and $G_{\rm RP}$, together with a universal HRD (the prior information on stellar models) can be used to improve distance estimates. Using the full spectra should ultimately improve this further.

Let $\mathbf{X} = (\varpi, G, G-G_{\rm BP}, G-G_{\rm RP}, G_{\rm BP}-G_{\rm RP})$ denote the set of measurements we have.
From these we want to determine the distance, but the measurements are also affected by the absolute magnitude, $M_V$, the intrinsic colour, $(V-I_C)$, as well as the interstellar extinction, $A_V$, so we will need to marginalise over these. We denote these parameters as $\mathbf{\Theta} = (r, A_V, M_V, (V-I_C))$. 
Bayes' theorem tells us that the posterior PDF over these parameters is
\begin{equation}
P(\mathbf{\Theta}|\mathbf{X}) = \frac{1}{Z}P(\mathbf{X}|\mathbf{\Theta})P(\mathbf{\Theta}),
\label{eq:post_mwhrd}
\end{equation}
where $Z$ is the normalisation constant, $P(\mathbf{X}|\mathbf{\Theta})$ is the likelihood to observe the data given the parameters, and $P(\mathbf{\Theta})$ is the prior probability of the parameters. 
The marginalised posterior over a single parameter of interest $\Theta_i$ is
\begin{equation}
P(\Theta_i|\mathbf{X}) = \int P(\mathbf{\Theta}|\mathbf{X}) \, d\Theta_1\ldots d\Theta_{i-1}d\Theta_{i+1}\ldots d\Theta_n \ .
\end{equation}
For simplicity, we consider the priors to be independent
\begin{equation}
P(\mathbf{\Theta}) = P(r)P(A_V)P((V-I_C),M_V) \ .
\end{equation}
$P((V-I_C),M_V)$ is the HRD. One could of course adopt more complex priors and additional parameters.

\subsection{The likelihood and forward model}

We will adopt a covariant Gaussian likelihood with mean $\mathbf{X}$ and covariance $\mathbf{\Sigma}$, i.e.\ \begin{equation}
P(\mathbf{X}|\mathbf{\Theta}) \propto \exp\left(-\frac{1}{2}\left[\mathbf{X}-\mathbf{f}(\mathbf{\Theta})\right]^T\mathbf{\Sigma}^{-1}\left[\mathbf{X}-\mathbf{f}(\mathbf{\Theta})\right]\right) \ .
\end{equation}
$\mathbf{f}(\mathbf{\Theta})$ is the forward model predicting the observations given the parameters:
\begin{equation}
\mathbf{f}(\mathbf{\Theta}) = 
\left(
\begin{array}{c}
1/r\\[0.5em]
M_V + 5\log_{10}r - 5 + A_V + p_0\left[(V-I_C)_{\rm obs}\right]\\[0.5em]
p_1\left[(V-I_C)_{\rm obs}\right]\\[0.5em]
p_2\left[(V-I_C)_{\rm obs}\right]\\[0.5em]
p_3\left[(V-I_C)_{\rm obs}\right]
\end{array}
\right) \ .
\end{equation}
$(V-I_C)_{\rm obs}$ is the observed colour, which is related to the intrinsic $(V-I_C)$ colour (parameter) by
\begin{equation}
(V-I_C)_{\rm obs} = (V-I_C) + A_V\left(1-\frac{A_I}{A_V}\right) \ ,
\end{equation}
where the ratio $A_I/A_V = 0.565$ was calculated using the \cite{fit99} extinction curve and is assumed to be constant throughout. The $p_i$ functions are polynomials which transform the $(V-I_C)$ colours into colours in the \textit{Gaia}{} photometric system, and are of the form
\begin{equation}
p_i[(V-I_C)] = c_0 + c_1(V-I_C) + c_2(V-I_C)^2 + c_3(V-I_C)^3 \ .
\end{equation}
The values of the coefficients are taken from \cite{jor10} and are tabulated in Table~\ref{tab:poly}.
\begin{table}
\caption{The coefficients of the transformation polynomials $p_i[(V-I_C)]$, from \cite{jor10}.}
\label{tab:poly}
\begin{center}
\begin{tabular}{clrrrr}
\hline\hline
\multicolumn{1}{c}{$i$} & \multicolumn{1}{c}{$p_i[(V-I_C)]$} & \multicolumn{1}{c}{$c_0$} &\multicolumn{1}{c}{$c_1$} &\multicolumn{1}{c}{$c_2$} &\multicolumn{1}{c}{$c_3$}\\
\hline
0 & $G-V$ & -0.0257 & -0.0924 & -0.1623 & 0.0090\\
1 & $G-G_{\rm BP}$ & 0.0387 & -0.4191 & -0.0736 & 0.0040\\
2 & $G-G_{\rm RP}$ & -0.0274 & 0.7870 & -0.1350 & 0.0082\\
3 & $G_{\rm BP}-G_{\rm RP}$ & -0.0660 & 1.2061 & -0.0614 & 0.0041\\
\hline
\end{tabular}
\end{center}
\end{table}

The covariance matrix $\mathbf{\Sigma}$ is constructed assuming that the parallax and photometric measurements are independent. The three photometric measurements are independent, but the colours formed from them are not. Thus the covariance matrix is
\begin{equation}
\mathbf{\Sigma} = 
\left(
\begin{array}{ccccc}
\sigma_\varpi^2 & 0 & 0 & 0 & 0\\[0.5em]
0 & \sigma_G^2 & \sigma_G^2 & \sigma_G^2 & 0\\[0.5em]
0 & \sigma_G^2 & \sigma_G^2+\sigma_{G_{\rm BP}}^2 & \sigma_G^2 & -\sigma_{G_{\rm BP}}^2\\[0.5em]
0 & \sigma_G^2 & \sigma_G^2 & \sigma_G^2+\sigma_{G_{\rm RP}}^2 & \sigma_{G_{\rm RP}}^2\\[0.5em]
0 & 0 & -\sigma_{G_{\rm BP}}^2 & \sigma_{G_{\rm RP}}^2 & \sigma_{G_{\rm BP}}^2+\sigma_{G_{\rm RP}}^2\\[0.5em]
\end{array}
\right) \ .
\end{equation}
To ensure that this matrix is not singular, we use regularisation for the inversion. We found that an offset of the order of 0.1 milimagnitudes in the photometric errors suffices, and has no relevant impact on the results.

\begin{figure}
\includegraphics[width=\hsize]{./plot_GUMS_KDE-IntrinsicHRD_2dHistogram_V-I_5sigma.pdf}
\caption{The HRD prior $P((V-I_C),M_V)$. The colour scale indicates the probability density and is normalised to unity. The smallest value is $1.25\times 10^{-9}$,
which is the value (post-normalization) of the offset added to avoid zero probability.}
\label{fig:hrd}
\end{figure}

\subsection{The priors}
For the prior on distance $P(r)$ we use again the exponentially decreasing space density{} prior discussed in Sect.~\ref{subsec:3priors}, using $L=1.35$~kpc, as well as the Milky Way density model $\rho_{\rm MW}(r,l,b)$ discussed in Sect.~\ref{subsec:mw}. The extinction prior $P(A_V)$ is a Gaussian centered on $A_V^{\rm map}(r,l,b)$ with standard deviation $\sigma_{A_V} = 0.1A_V^{\rm map}$, where $A_V^{\rm map}(r,l,b)$ is the extinction calculated from the \cite{dri01} extinction map shown in Fig.~\ref{fig:extMap}.

The HRD prior $P((V-I_C),M_V)$ is constructed from the same GUMS data used to construct the luminosity function described in Sect.~\ref{subsec:mw}. Again we combine the CMDs of different components of the Milky Way into a single universal CMD (or HRD), and smooth it using two-dimensional kernel density estimation (a Gaussian kernel with widths of 0.05\,mag in both directions and no covariance is used). The number of stars per volume element is again used as the weighting to calculate the density of points in the HRD on a $2000\times2000$ grid. A small positive constant is added to ensure no region of the HRD has exactly zero probability density. The value of this is set such that its ratio to the total probability (original plus offset density) is $2.87\times 10^{-7}$ (this number corresponds to 5$\sigma$ probability in a Gaussian distribution). 
This total probability is then normalised to unity. The resulting HRD prior is shown in Fig.~\ref{fig:hrd}. It it far from perfect, but it is adequate for demonstration purposes.

\begin{figure}
\includegraphics[width=\hsize]{./plot_GUMS_09_photometricError_small.pdf}
\caption{The distribution of the end-of-mission photometric errors $\sigma_{G}$ (left), $\sigma_{G_{\rm BP}}$ (middle), and $\sigma_{G_{\rm RP}}$ (right) as a function of true $G$ magnitude, for all stellar objects in the GUMS catalogue. The colour shows the probability density per unit magnitude per unit log photometric errors.}
\label{fig:photError}
\end{figure}
\subsection{The photometric error model}
The error models for the photometric measurements are the parameterised form\footnote{\url{www.cosmos.esa.int/web/gaia/science-performance}} of the formulation outlined in \cite{jdb022} and \cite{jor10}. This includes the scattered light levels measured post-launch \citep{deb14}. Using the same notation as for our parallax noise model (Eq.~\ref{eq:sigmaVarpi}), the standard errors in the $G$-band magnitude is
\begin{multline}
\sigma_G\;[{\rm mag}] = m\left[\frac{1}{n_{\rm tr}p_{{\rm det},G}(G)}\left(\sigma_{\rm cal}^2 \right.\right.\\ 
\left.\vphantom{\frac{1}{n_{\rm tr}p_{{\rm det},G}(G)}}\left.\vphantom{\sigma_{\rm cal}^2} + 10^{-6}\left[0.0001985 + 1.8633z + 0.04895z^2\right]\right)\right]^{1/2},
\label{eq:sigmaG}
\end{multline}
where
\begin{equation}
z = \max\left(10^{0.4(12 - 15)}, 10^{0.4(G-15)}\right),
\end{equation}
$m = 1.2$ is a contingency margin to account for unknown sources of error, and $\sigma_{\rm cal} = 30$\,mmag \citep{jdb022, jor06, cj047} is the adopted total photocalibration error.

\begin{table}
\caption{The coefficients $c_{{\rm XP},ij}$ used in Eqs.~\ref{eq:sigmaGbprp}--\ref{eq:coeffXp} to calculate the BP/RP photometric errors $\sigma_{G_{\rm XP}}$.}
\begin{tabular}{clrrrr}
\hline\hline
\multicolumn{1}{c}{XP} & \multicolumn{1}{c}{$c_{{\rm XP},i}$} & \multicolumn{1}{c}{$c_{{\rm XP},i0}$} &\multicolumn{1}{c}{$c_{{\rm XP},i1}$} &\multicolumn{1}{c}{$c_{{\rm XP},i2}$} &\multicolumn{1}{c}{$c_{{\rm XP},i3}$}\\
\hline
\multirow{3}{*}{BP} & $c_{{\rm BP},0}$ & -1.866968 & -0.205807 & 0.060769 & 0.000262\\
& $c_{{\rm BP},1}$ & 1.465592 & 0.195768 & 0.018878 & -0.000400\\
& $c_{{\rm BP},2}$ & 1.043270 & 0.355123 & 0.044390 & -0.000562\\
\hline
\multirow{3}{*}{RP} & $c_{{\rm RP},0}$ & -3.042268 & -0.091569 & 0.027352 & -0.001923\\
& $c_{{\rm RP},1}$ & 1.783906 & -0.318499 & 0.057112 & -0.003803\\
& $c_{{\rm RP},2}$ & 1.615927 & -0.636628 & 0.114126 & -0.007597\\
\hline
\end{tabular}
\label{tab:cxp}
\end{table}
The BP/RP photometric errors are
\begin{multline}
\sigma_{G_{\rm XP}}\;[{\rm mag}] = \left[\frac{1}{n_{\rm tr}p_{{\rm det},XP}(G)}\right.\times\\
 \left.\vphantom{\frac{1}{n_{\rm tr}p_{{\rm det},XP}(G)}}
 \left(1.44\sigma_{\rm cal}^2 + 10^{-6}\sum_{i=0}^{2}10^{c_{{\rm XP},i}}z^i\right)\right]^{1/2}
\label{eq:sigmaGbprp}
\end{multline}
where
\begin{equation}
z = \max\left(10^{0.4(11 - 15)}, 10^{0.4(G-15)}\right),
\end{equation}
the coefficients $c_{{\rm XP},i}$ are polynomials in $(V-I_C)$:
\begin{equation}
c_{{\rm XP},i} = \sum_{J=0}^{3}c_{{\rm XP},ij}(V-I_C)^j,
\label{eq:coeffXp}
\end{equation}
in which the values of the coefficients $c_{{\rm XP},ij}$ are tabulated in Table~\ref{tab:cxp}, and the BP/RP detection probability $p_{\rm XP}(G)$ is taken from Table~10 in \cite{jor10}. The BP/RP parameterisation in Eqs.~\ref{eq:sigmaGbprp}--\ref{eq:coeffXp} already includes the contingency margin of $m=1.2$.

The distributions of end-of-mission $\sigma_G$, $\sigma_{G_{\rm BP}}$, and $\sigma_{G_{\rm RP}}$ as a function of the $G$-band magnitude are shown in Fig.~\ref{fig:photError}. BP photometric errors $\sigma_{G_{G_{\rm BP}}}$ shows larger errors than the other two bands because of the effects of extinction. We can see in Fig.~\ref{fig:photError} that there is an additional lower density locus in the distribution of $G_{\rm BP}$, which is due to areas with heavy extinction in the GUMS catalogue.

To simulate the photometric observations from GUMS, we draw the $G$ and $G_{\rm XP}$ magnitudes 
as we did the parallaxes, by drawing at random from their likelihoods (which are independent 1D Gaussians with mean equal to the true magnitude and standard deviation given by the noise model).

\subsection{Computing the posterior}

Since the posterior is now multivariate, we now use the affine-invariant Markov Chain Monte Carlo (MCMC) ensemble sampler \citep{goo10} to sample the multivariate posterior, from which the one-dimensional marginalised posterior $P(r|\mathbf{X})$ is easy to obtain via density estimation over these samples. For each source, samples are taken in a maximum of 30\,000 steps, but for every 10\,000 steps the autocorrelation time is calculated and the chain is stopped when convergence is achieved.

The extinction $A_V^{\rm map}(r,l,b)$ is calculated using a lookup table that tabulates the extinction in distance and \texttt{healpix} cell, with $n_{\rm side} = 64$. For a source with a given $(l,b)$ we find the cell closest to that direction and interpolate $A_V$ in distance. This speeds up the calculation in two ways: We do not have to calculate the extinction using the full set of equations and parameters that describe the extinction model, and we can reuse the extinction curve in a given cell for all stars in the cell (there are far fewer cells than stars).

\begin{figure}
\includegraphics[width=\hsize]{./plot_GUMS_comparison_APTrueVSAPMo_05_priorsWithMWHRDMCMC_logBin.pdf}
\caption{Comparison of the mode parameters, from top to bottom: Distance $r\mode$, extinction $A_{V,Mo}$, absolute magnitude $M_{V,Mo}$, and colour $(V-I_C)_{Mo}$ for all stars in the GUMS catalogue, with their corresponding true values, inferred using the exponentially decreasing space density{} prior (left column) and the Milky Way{} (right column) prior as the distance prior. The black lines in each panel indicate a perfect match between the estimated and true value.}
\label{fig:comparisons_mwhrd_mo}
\end{figure}
\begin{figure*}
\centering
\includegraphics{plot_GUMS_summaryPerformance_mode_longRangeF_03_3panels.pdf}
\caption{Performance of the mode distance $r\mode$ using the exponentially decreasing space density\,+\,HRD\,+\,phot{} prior (thick teal and cyan lines) and the Milky Way\,+\,HRD\,+\,phot{} prior (thick violet and pink lines), shown in terms of the bias $\overline{x}_{Mo}$ (top rows), the RMS (middle rows), and the standard deviation $\sigma_{x,\mode}$ (bottom rows), as a function of $f\true$ (left columns) and $f_{\rm obs}$ (right columns). 
For both distance priors, the lines in lighter shades (viz. cyan and pink) are for all stars, whereas the lines in darker shades (teal and violet) are just for stars with photometric errors less than 50\,mmag.
The performance of the corresponding parallax-only distance estimates (thin blue and magenta lines) previously shown in Fig.~\ref{fig:summary_performance_ftrue_mode} and Fig.~\ref{fig:summary_performance_fobs_mode} are included for comparison. Note that the vertical axes have different (smaller) ranges than in Fig.~\ref{fig:summary_performance_ftrue_mode} and Fig.~\ref{fig:summary_performance_fobs_mode}.}
\label{fig:summary_performance_ftrue_mode_mwhrd}
\end{figure*}
\begin{figure*}
\includegraphics[width=\hsize]{plot_GUMS_summaryPerformanceFTrue_mode_longRangeF_directional_RMS_03_combined.pdf}
\caption{As Fig.~\ref{fig:summary_performance_ftrue_mode_dir}, but for the exponentially decreasing space density\,+\,HRD\,+\,phot{} prior and the Milky Way\,+\,HRD\,+\,phot{} prior. The results of the corresponding parallax-only distance estimates are also included for comparison. Note that the vertical axes have different (smaller) ranges than in Fig.~\ref{fig:summary_performance_ftrue_mode_dir}.}
\label{fig:summary_performance_ftrue_mode_mwhrd_dir}
\end{figure*}
\begin{figure}
\includegraphics[width=\hsize]{plot_GUMS_formalErrors_ZHistogram_rMo_03_combined.pdf}
\caption{As Fig.~\ref{fig:zhist}, but for the exponentially decreasing space density\,+\,HRD\,+\,phot{} prior and the Milky Way\,+\,HRD\,+\,phot{} prior, for all stars and with the photometric error cutoff applied. The results for the corresponding parallax-only inferences are shown for comparison. The black dashed line indicates the standard normal distribution. All distributions are scaled to have their maxima at one. Note the different range on the $z$-axis compared to Fig.~\ref{fig:zhist}.}
\label{fig:zhist2}
\end{figure}

\subsection{Results of the joint astrometric and photometric posterior}

Having computed the posterior we find the mode of the parameter vector, $\mathbf{\Theta}_{Mo}$, for each star, by finding the sample in the posterior chain that gives the highest posterior probability.
Fig.~\ref{fig:comparisons_mwhrd_mo} compares this with the true parameters, $\mathbf{\Theta}_{\rm true}$, for both distance priors. We see that the distance inference---our main goal---is much improved compared to the distance estimate using just the parallax and (the same) distance priors shown in Fig.~\ref{fig:comparisons_mo}. Comparing the two distance priors (left and right column in Fig.~\ref{fig:comparisons_mwhrd_mo}), we see no significant difference, which indicates that any poor distance estimates arising from poorly measured parallaxes can be improved by the inclusion of photometric data.

Despite the improvement, we see that for $r\true\gtrsim 2$\,kpc, the inferred distances are underestimated for some stars. This is primarily because we are unable to distinguish dwarfs from giants with these colours. Looking at the inference for the absolute magnitudes, we see that the region with $M_{V,{\rm true}}$ between about -3 and 5 is where the absolute magnitudes are slightly overestimated, i.e.\ they are inferred to be fainter than the true values. The model is therefore placing these stars closer than they really are.

The performance of the distance estimation as a function of $f\true$ and $f_{\rm obs}$ averaged over all directions is shown in Fig.~\ref{fig:summary_performance_ftrue_mode_mwhrd}.
The primary comparison here is between the line labelled ``exponentially decreasing space density\,+\,HRD\,+\,phot{}'' and line labelled ``exponentially decreasing space density{}'', as well as between the line labelled ``Milky Way\,+\,HRD\,+\,phot{}'' and the line labelled ``Milky Way{}''. Looking at the plots against $f_{\rm obs}$ first (right column), we see the significant improvement in performance in both the bias and standard deviation when adding the photometry (and the HRD prior). Indeed, when using the Milky Way{} model as the distance prior, we see that adding the photometry almost entirely eliminates the bias for small positive $f_{\rm obs}$ and stabilises it at $\overline{x}_{Mo}\sim 0.08$ for $f_{\rm obs}\gtrsim 1$. Using the exponentially decreasing space density{} prior eliminates the bias even further, with $\overline{x}_{Mo}\sim -0.03$ for the same range of $f_{\rm obs}$. The negative value indicates that the distances to most stars are still slightly underestimated, due to the constraint imposed by the distance prior.
The inclusion of photometry also significantly reduces the standard deviation in the distance estimates, although interestingly, with the Milky Way{} model as the distance prior, it does not decrease the standard deviation more than what we achieved using the exponentially decreasing space density{} prior with just the parallax data.

Nonetheless, the use of parallaxes and photometry together produces smaller RMS residuals than using the parallax alone (with any of the priors). We can see this in the right column and middle row of Fig.~\ref{fig:summary_performance_ftrue_mode_mwhrd}. Using the Milky Way{} prior, the RMS reaches a plateau of around $\sim$0.45 with a bias of $\sim$0.08 once $f_{\rm obs}$ rises above $\sim$1. For negative parallaxes, the bias and the RMS are stable at respectively $\sim$0.08 and $\sim$0.42. In this regime, the inclusion of photometric data makes this prior superior to other parallax-only priors. Restricting the analysis to stars with small photometric uncertainties (less than 50\,mmag in all bands) only reduces the bias by a negligible amount, but improves the RMS errors and standard deviation by $\sim$0.05 for both positive and negative $f_{\rm obs}$. Further improvements can be achieved when the exponentially decreasing space density{} prior is instead used as the distance prior. For positive parallaxes, the bias is only about $\sim$-0.03, and the RMS is about $\sim$0.42. Considering only stars with small photometric uncertainties (the same cut as before) improves these further by about the same amount as using the Milky Way{} prior. Having said this, it should be recalled that the majority of stars with positive parallaxes and $f_{\rm obs}\lesssim 3$ are more likely to be disc stars. This is why the exponentially decreasing space density\,+\,HRD\,+\,phot{} prior performs better than the Milky Way\,+\,HRD\,+\,phot{} prior. 

In terms of $f\true$ (left column of Fig.~\ref{fig:summary_performance_ftrue_mode_mwhrd}), we see that the Milky Way\,+\,HRD\,+\,phot{} prior is superior to the parallax-only Milky Way{} prior, although over a range of $f\true$ between 0.1 and 1.8 the parallax-only exponentially decreasing space density{} prior is actually better in RMS. Using only stars with small photometric uncertainties improves the distance estimates, with the RMS and standard deviation are always below $\sim$0.45 for all $f\true$. If we look at the distribution of the photometric errors in Fig.~\ref{fig:photError}, we see that the 50\,mmag cutoff affects mostly red-colored stars (whether this be intrinsic or due to reddening). We see that the exponentially decreasing space density\,+\,HRD\,+\,phot{} prior improves the estimation further, with the RMS always below $\sim$0.4 even when all stars are included. Taking only stars with small photometric uncertainties reduces the RMS further by $\sim$0.05. For poorly measured stars with $f\true\gtrsim 1.2$, however, the Milky Way\,+\,HRD\,+\,phot{} prior performs better with the RMS stabilising at $\sim$0.3 for large $f\true$. We shall see later that this is due to the poor performance of the exponentially decreasing space density\,+\,HRD\,+\,phot{} prior for directions toward the galactic centre.

The directional performance (in terms of the RMS) of the mode estimator for the exponentially decreasing space density\,+\,HRD\,+\,phot{} and the Milky Way\,+\,HRD\,+\,phot{} prior is shown in Fig.~\ref{fig:summary_performance_ftrue_mode_mwhrd_dir}. Here we see that the performance is relatively stable everywhere, both priors performs more-or-less equally well in directions away from the galactic plane, with the exponentially decreasing space density\,+\,HRD\,+\,phot{} performing $\sim$0.05 better for large $f\true$. Both are especially accurate in the galactic plane towards the direction of the anticentre (at $l=0^\circ$ and $b=135^\circ, 180^\circ$), where the RMS stabilise at $\sim$0.2 for large $f\true$ (in other directions the RMS stabilise at $\sim$0.25--0.4 for large $f\true$). In some directions, however, some ``spikes'' can be seen in the performance curve of the Milky Way\,+\,HRD\,+\,phot{} prior. These corresponds to some stars which distances are so poorly estimated they disrupt the overall behaviour of the curve. The distances are poorly estimated again because of the inability to distinguish dwarfs from giants. These spikes tend to appear at high latitudes due to the paucity of stars in these directions. At the galactic plane, the stars are more numerous, the overall behaviour is relatively undisturbed by these stars. Making cuts in photometric errors as before can remove most of these stars and improves the overall performance, suggesting that these poor performances are to some degree due to poor photometric measurements. This shows that a combination of bad photometric and parallax measurements will not help in accurately estimating distances. Bad parallax measurements should be complemented by good photometric measurements and vice versa.

These directions aside, the galactic centre (bottom-right panel of Fig.~\ref{fig:summary_performance_ftrue_mode_mwhrd_dir}) remains a problematic direction. For the Milky Way\,+\,HRD\,+\,phot{} prior, we see here that the RMS error in this direction can reach up to 1.64 even for $f\true\lesssim 0.5$. The performance improves as $f\true$ increase, but it is hardly better than the parallax-only Milky Way{} prior. Considering only stars with small photometric uncertainties improves the situation considerably, although it is not better than the exponentially decreasing space density{} prior even for $f\true\lesssim 1$. For the same range of $f\true$, the exponentially decreasing space density\,+\,HRD\,+\,phot{} prior performs best for this direction, but only up to $f\true\sim 1.2$. For larger $f\true$, the Milky Way\,+\,HRD\,+\,phot{} prior (with cuts in photometric uncertainties) performs better. The heavy extinction in this direction seem to be main cause of the poor performance of the Milky Way\,+\,HRD\,+\,phot{} prior, but for the exponentially decreasing space density\,+\,HRD\,+\,phot{} prior the poor performance is due to the mismatch between the distance prior and true distribution. Clearly $L=1.35$\,kpc is not a good scale length for directions toward the centre. Since the Milky Way{} model takes the galactic bulge into account (see Fig~\ref{fig:prior_mw} for the same direction), it provides a better distance constraint. Considering the performance in this direction, we can see that the averaged performance shown in Fig.~\ref{fig:summary_performance_ftrue_mode_mwhrd} is dominated by stars towards the galactic centre. 

Finally, we compare the formal uncertainties (Sect.~\ref{subsec:formalErrors}) in Fig.~\ref{fig:zhist2}. The distributions of the standardised residuals $z$ for the two pairs are similar, an indication that the uncertainties in the distance estimation depends largely on photometric data. The $z$-distribution is asymmetric, although for positive $z$ the distribution is very close to a standard normal distribution (shown with a black dashed line). For the underestimated distances (negative $z$), the width is less than that of a standard normal distribution, suggesting that the uncertainties are more overestimated than the uncertainties for the overestimated distances (positive $z$). Overall, the HWHM of the $z$-distribution is 0.9, which suggests that in general the uncertainties might be overestimated, altough not as much as the parallax-only Milky Way{} prior (which is 0.84). The HWHM of the HRD + phot priors seems to suggest that the inclusion of photometric data can produce a smaller credible interval on distance (hence a larger value of the HWHM), but the photometric errors also add additional uncertainties thus keeping it smaller than the HWHM of a standard Gaussian.

\section{Summary and conclusions}
We have estimated distances from simulated \textit{Gaia}{} parallaxes using probability-based inference. We examined four different priors on distances: the three isotropic priors---uniform distance{}, uniform space density{}, exponentially decreasing space density{}---and the more complex Milky Way{} prior, which models the Galaxy and its observability in a magnitude-limited survey. We further considered the inclusion of photometric measurements (plus an HRD prior) with the exponentially decreasing space density{} prior and the Milky Way{} prior. The simulated \textit{Gaia}{} data were taken from the \textit{Gaia}{} Universe Model Snapshot (GUMS) catalogue. Note that the Milky Way{} prior was intentionally not taken from this. We investigated the use of both the median and the mode of the posterior as distance estimators. The median was found to be most appropriate for the uniform space density{} prior, but the mode performed better for the other priors.

For stars with positive parallaxes and observed fractional parallax errors $f_{\rm obs}$ less than about 0.1, the distance estimates were largely independent of the choice of prior. This will only apply to about $\sim$10\% of the \textit{Gaia}{} catalogue, however. Taking this approach would not only limit us to a small fraction of hard-won data, it would also lead to truncation biases when using such a sample for astrophysical analyses.

We found that the uniform space density{} prior with an upper limit of $r_{\rm lim}=100$\,kpc performs poorly. This is partly because this is an unrealistic prior when applied to a large volume of the Galaxy, since it is obvious that within a radius of $r=100$\,kpc the volume density of the Galaxy is not constant (the mode should not be used with this prior because once $f_{\rm obs}$ is larger than about 0.35, the mode is at $r_{\rm lim}$). Yet truncating this prior at smaller distances produces other, no less dramatic, problems.

Viewed as a function of the expected (true) fractional parallax error, $f\true$, the uniform distance{} prior also does very poorly, showing a steep rise in bias and standard deviation above $f\true=0.1$ (already over 100\% at $f\true=0.3$).
But if we view it as a function of the observed fractional parallax error, $f_{\rm obs}$, we see much better results. This is because plotting against $f_{\rm obs}$ automatically ``cleans'' the appearance by assigning poor data (very small or negative parallaxes) to large or negative values of $f_{\rm obs}$.
Nonetheless, over moderate positive values of $f_{\rm obs}$ this prior works reasonably well, better than the Milky Way{} prior and not much worse than the exponentially decreasing space density{} prior with $L=1.35$\,kpc. However, a closer inspection shows that most of its good performance comes from a narrow region towards the galactic centre where there are a lot of stars (see Fig.~\ref{fig:summary_performance_fobs_mode_dir}), which therefore dominates the average performance. The performance at other low latitude directions is considerably worse than the Milky Way{} and exponentially decreasing space density{} priors, and at high latitudes it is no better than the exponentially decreasing space density{} (and significantly worse once $f_{\rm obs}\gtrsim 1$).

The exponentially decreasing space density{} prior performs well, altough fine-tuning is required to find the appropriate scale length $L$ (the only parameter in the prior). The RMS of the scaled residual in distance (Eq.~\ref{eq:x}), averaged over all directions in the galaxy, is roughly equal to $f\true$ for $f\true\lesssim 0.4$, but does not increase further as $f\true$ grows to 1. In terms of these RMS errors, we found that $L=1.35$\,kpc gives the best performance. The mode of this prior is $2L = 2.7$\,kpc, and this is on the order of the scale length of the stellar disc of the Milky Way, where most stars are located. Thus it can accurately infer the distance of a large majority of stars in the catalogue. We see, however, that the RMS increases with $f\true$, because poor parallaxes ``force'' the mode of the posterior to be close to the mode of the prior at $2L$. This is far from the true distances for many of these stars, as they are more likely to be distant stars.

The Milky Way{} prior performs well in directions away from the galactic centre. The poor performance toward the galactic centre is due a prior mismatch problem, specifically due to a different adopted bulge-to-disc mass ratio for our prior. This can be seen from a comparison between its distance dependence and that in the ``true'' Galaxy: We see in Fig.~\ref{fig:gums_rdist_directional} significant discrepancies at directions close to the galactic centre. Again, as in the case of the uniform distance{} prior, the performance in this direction dominates the direction-averaged performance. The prior could be improved by adopting a different mass density model and by using a different luminosity function $\phi(M_G)$ for each component of the Milky Way, rather than a single universal luminosity function. However, such discrepancies between the prior and the underlying true distribution we are trying to infer will inevitably occur in any real application. Our results show that although constructing a detailed model of the Milky Way is a good idea in principle, in practice it can go wrong due to the strong assumptions involved and subsequent likelihood of a mismatch.

Rather than creating a detailed prior, a better way to improve distance estimation is to include other distance-sensitive measurements, where possible. We have shown that photometric data can help for stars with poor parallax measurements. This requires us to make further assumptions, namely about the intrinsic nature of the sources. To make use of photometry we (must) adopt a Hertzsprung--Russell diagram and an extinction map. We find that this improves the averaged performance significantly, with an RMS error which is always less than about 0.5 for $f\true\lesssim 3$. It is, however, still difficult to estimate distances for stars toward the galactic centre, due to the heavy extinction in this direction. At higher latitudes the performance is more or less stable.

While adding photometry helps, \textit{Gaia}{} broad band photometry on its own does not allow us to distinguish giants from dwarfs, so this remains a limitation in the distance accuracy even when including photometric data. One could of course use more sophisticated priors or more data such as spectroscopy (e.g. \citet{sch14, san16, liu12}) to help estimate absolute luminosity, but one must be aware that distance estimates which do not use just the parallax necessarily require a different set of assumptions, such as the physical properties source. The stronger the prior, the greater the risk of mismatch problems. Increased use of priors may reduce the variance in the distance estimations, but they can also increase the bias.

\acknowledgements
We thank Jos de Bruijne for the discussion on \textit{Gaia}{} astrometric errors, Berry Holl for providing us with the transit data in Fig.~\ref{fig:ntr}, the DPAC/CU2 developers of \texttt{GaiaSimu} for the code to calculate the extinction map in Fig.~\ref{fig:extMap}, Joe Hennawi and the Max Planck Computing and Data Facility (MPCDF) in Garching for the computational resources, and the IT group of the Max Planck Institute for Astronomy for their technical assistance. We thank Anthony G.A. Brown, Morgan Fouesneau, and Jan Rybizki for useful comments and discussions. We also thank the anonymous referee for the comments and suggestions that improve the analysis of this paper.

This research made use of NASA's Astrophysics Data System; \texttt{Astropy}, a community-developed core Python package for astronomy \citep{2013A&A...558A..33A}; the \texttt{IPython} package \citep{PER-GRA:2007}; \texttt{matplotlib}, a Python library for publication quality graphics \citep{Hunter:2007}; and \texttt{SciPy} \citep{jones_scipy_2001}. 

