  The $^{16}$O(p,$\gamma$)$^{17}$F reaction is the slowest hydrogen-burning
process in the CNO mass region. Its thermonuclear rate sensitively impacts
predictions of oxygen isotopic ratios in a number of astrophysical sites,
including AGB stars. The reaction has been measured several times at low
bombarding energies using a variety of techniques. The most recent evaluated
experimental rates have a reported uncertainty of about 7.5\% below $1$~GK.
However, the previous rate estimate represents a best guess only and was not
based on rigorous statistical methods. We apply a Bayesian model to fit all
reliable $^{16}$O(p,$\gamma$)$^{17}$F cross section data, and take into account
independent contributions of statistical and systematic uncertainties. The
nuclear reaction model employed is a single-particle potential model involving
a Woods-Saxon potential for generating the radial bound state wave function.
The model has three physical parameters, the radius and diffuseness of the
Woods-Saxon potential, and the asymptotic normalization coefficients (ANCs) of
the final bound state in $^{17}$F. We find that performing the Bayesian $S$
factor fit using ANCs as scaling parameters has a distinct advantage over
adopting spectroscopic factors instead. Based on these results, we present the
first statistically rigorous estimation of experimental
$^{16}$O(p,$\gamma$)$^{17}$F reaction rates, with uncertainties ($\pm 4.2$\%)
of about half the previously reported values.
