  The amplitude of the ionizing background that pervades the intergalactic
medium (IGM) at the end of the epoch of reionization provides a valuable
constraint on the emissivity of the sources which reionized the Universe. While
measurements of the ionizing background at lower redshifts rely on a
simulation-calibrated mapping between the photoionization rate and the mean
transmission of the Ly$\alpha$ forest, at $z\gtrsim6$ the IGM becomes
increasingly opaque, and transmission arises solely in narrow spikes separated
by saturated Gunn-Peterson troughs. In this regime, the traditional approach of
measuring the average transmission over large $\sim 50$ Mpc$/h$ regions is less
sensitive and sub-optimal. Additionally, the five times smaller oscillator
strength of the Ly$\beta$ transition implies the Ly$\beta$ forest is
considerably more transparent at $z\gtrsim6$, even in the presence of
contamination by foreground $z\sim 5$ Ly$\alpha$ forest absorption. In this
work we present a novel statistical approach to analyze the joint distribution
of transmission spikes in the co-spatial $z\sim 6$ Ly$\alpha$ and Ly$\beta$
forests. Our method relies on Approximate Bayesian Computation (ABC), which
circumvents the necessity of computing the intractable likelihood function
describing the highly correlated Ly$\alpha$ and Ly$\beta$ transmission. We
apply ABC to mock data generated from a large-volume hydrodynamical simulation
combined with a state-of-the-art model of ionizing background fluctuations in
the post-reionization IGM, and show that it is sensitive to higher IGM neutral
hydrogen fractions than previous techniques. As a proof of concept, we apply
this methodology to a real spectrum of a $z=6.54$ quasar and measure the
ionizing background from $5.4\leq z \leq 6.4$ along this sightline with
$\sim0.2$ dex statistical uncertainties.
