  Observations of thermonuclear X-ray bursts from accreting neutron stars (NSs)
in low-mass X-ray binary systems can be used to constrain NS masses and radii.
Most previous work of this type has set these constraints using Planck function
fits as a proxy: both the models and the data are fit with diluted blackbody
functions to yield normalizations and temperatures which are then compared
against each other. Here, for the first time, we fit atmosphere models of X-ray
bursting NSs directly to the observed spectra. We present a hierarchical
Bayesian fitting framework that uses state-of-the-art X-ray bursting NS
atmosphere models with realistic opacities and relativistic exact Compton
scattering kernels as a model for the surface emission. We test our approach
against synthetic data, and find that for data that are well-described by our
model we can obtain robust radius, mass, distance, and composition
measurements. We then apply our technique to Rossi X-ray Timing Explorer
observations of five hard-state X-ray bursts from 4U 1702-429. Our joint fit to
all five bursts shows that the theoretical atmosphere models describe the data
well but there are still some unmodeled features in the spectrum corresponding
to a relative error of 1-5% of the energy flux. After marginalizing over this
intrinsic scatter, we find that at 68% credibility the circumferential radius
of the NS in 4U 1702-429 is R = 12.4+-0.4 km, the gravitational mass is
M=1.9+-0.3 Msun, the distance is 5.1 < D/kpc < 6.2, and the hydrogen mass
fraction is X < 0.09.
