Modern surveys of gravitational microlensing events have progressed to
detecting thousands per year. Surveys are capable of probing Galactic
structure, stellar evolution, lens populations, black hole physics, and the
nature of dark matter. One of the key avenues for doing this is studying the
microlensing Einstein radius crossing time distribution (tE). However,
systematics in individual light curves as well as over-simplistic modeling can
lead to biased results. To address this, we developed a model to simultaneously
handle the microlensing parallax due to Earth's motion, systematic instrumental
effects, and unlensed stellar variability with a Gaussian Process model. We
used light curves for nearly 10,000 OGLE-III and IV Milky Way bulge
microlensing events and fit each with our model. We also developed a forward
model approach to infer the timescale distribution by forward modeling from the
data rather than using point estimates from individual events. We find that
modeling the variability in the baseline removes a source of significant bias
in individual events, and previous analyses over-estimated the number of long
timescale (tE>100 days) events due to their over simplistic models ignoring
parallax effects and stellar variability. We use our fits to identify hundreds
of events that are likely black holes.Comment: Submitted version, in review, 33 pages, 18 figures, MCMC posterior
samples available by publisher after acceptanc