We develop a novel Bayesian methodology aimed at reliably and precisely inferring the distribution of dark matter within the Milky Way using rotation curve data. We identify a subset of the available rotation curve tracers that are mutually consistent with each other, thus eliminating data sets that might suffer from systematic bias. We investigate different models for the mass distribution of the luminous (baryonic) component that bracket the range of likely morphologies. We demonstrate the statistical performance of our method on simulated data in terms of coverage, fractional distance, and mean squared error. Applying it to Milky Way data we measure the local dark matter density at the solar circle ρ0 to be ρ0 = 0.43 ± 0.02(stat) ± 0.01(sys) GeV/cm3, with an accuracy ~ 6%. This result is robust to the assumed baryonic morphology. The scale radius and inner slope of the dark matter profile are degenerate and cannot be individually determined with high accuracy. We show that these results are robust to several possible residual systematic errors in the rotation curve data.