  We present a Bayesian method to constrain the masses and radii of neutron
stars (NSs) using the information encoded in the X-ray pulse profiles of
accreting millisecond pulsars. We model the shape of the pulses using "oblate
Schwarzschild" approximation, which takes into account the deformed shape of
the star together with the special and general relativistic corrections to the
photon trajectories and angles. The spectrum of the radiation is obtained from
an empirical model of Comptonization in a hot slab in which a fraction of seed
blackbody photons is scattered into a power-law component. By using an
affine-invariant Markov chain Monte Carlo ensemble sampling method, we obtain
posterior probability distributions for the different model parameters,
especially for the mass and the radius. To test the robustness of our method,
we first analyzed self-generated synthetic data with known model parameters.
Similar analysis was then applied for the observations of SAX J1808.4-3658 by
the Rossi X-ray Timing Explorer (RXTE). The results show that our method can
reproduce the model parameters of the synthetic data, and that accurate
constraints for the radius can be obtained using the RXTE pulse profile
observations if the mass is a priori known. For a mass in the range 1.5-1.8
Msun, the radius of the NS in SAX J1808.4-3658 is constrained between 9 and 13
km. If the mass is accurately known, the radius can be determined with an
accuracy of 5% (68% credibility). For example, for the mass of 1.7 Msun the
equatorial radius is Req = 11.9+0.5 -0.4 km. Finally, we show that further
improvements can be obtained when the X-ray polarization data from the Imaging
X-ray Polarimeter Explorer will become available.
