  We introduce a new approach for analyzing the IGM damping wings imprinted on
the proximity zones of quasars in the epoch of reionization (EoR). Whereas past
work has typically forgone the additional constraining power afforded by the
blue side continuum ($1216\,\r{A} \lesssim \lambda \lesssim 1280\,\r{A}$)
and/or opted not to model the large correlated IGM transmission fluctuations in
the proximity zone ($\lambda \lesssim 1216\,\r{A}$), we construct a generative
probabilistic model for the entire spectrum accounting for all sources of error
- the stochasticity induced by patchy reionization, the impact of the quasar's
ionizing radiation on the IGM, the unknown intrinsic spectrum of the quasar,
and spectral noise. This principled Bayesian method allows us to marginalize
out nuisance parameters associated with the quasar's radiation and its unknown
intrinsic spectrum to precisely measure the IGM neutral fraction, $\langle
x_{\rm HI}\rangle$. A key element of our analysis is the use of dimensionality
reduction (DR) to describe the intrinsic quasar spectrum via a small number of
nuisance parameters. Using a large sample of $15,559$ SDSS/BOSS quasars at $z
\gtrsim 2.15$ we trained and quantified the performance of six distinct DR
methods, and find that a six parameter PCA model (five coefficients plus a
normalization) performs best, with complex machine learning approaches
providing no advantage. By conducting statistical inference on 100 realistic
mock EoR quasar spectra, we demonstrate the reliability of the credibility
contours that we obtain on $\langle x_{\rm HI}\rangle$ and the quasar lifetime,
$t_{\rm Q}$. The new method introduced here will transform IGM damping wings
into a precision probe of reionization, on the same solid methodological and
statistical footing as other precision cosmological measurements.
