© 2016. The American Astronomical Society. All rights reserved. A principal scientific goal of the Gemini Planet Imager (GPI) is obtaining milliarcsecond astrometry to constrain exoplanet orbits. However, astrometry of directly imaged exoplanets is subject to biases, systematic errors, and speckle noise. Here, we describe an analytical procedure to forward model the signal of an exoplanet that accounts for both the observing strategy (angular and spectral differential imaging) and the data reduction method (Karhunen-Loève Image Projection algorithm). We use this forward model to measure the position of an exoplanet in a Bayesian framework employing Gaussian processes and Markov-chain Monte Carlo to account for correlated noise. In the case of GPI data on β Pic b, this technique, which we call Bayesian KLIP-FM Astrometry (BKA), outperforms previous techniques and yields 1σ errors at or below the one milliarcsecond level. We validate BKA by fitting a Keplerian orbit to 12 GPI observations along with previous astrometry from other instruments. The statistical properties of the residuals confirm that BKA is accurate and correctly estimates astrometric errors. Our constraints on the orbit of β Pic b firmly rule out the possibility of a transit of the planet at 10-σ significance. However, we confirm that the Hill sphere of β Pic b will transit, giving us a rare chance to probe the circumplanetary environment of a young, evolving exoplanet. We provide an ephemeris for photometric monitoring of the Hill sphere transit event, which will begin at the start of April in 2017 and finish at the end of January in 2018.