  A commonly used measure to summarize the nature of a photon spectrum is the
so-called Hardness Ratio, which compares the number of counts observed in
different passbands. The hardness ratio is especially useful to distinguish
between and categorize weak sources as a proxy for detailed spectral fitting.
However, in this regime classical methods of error propagation fail, and the
estimates of spectral hardness become unreliable. Here we develop a rigorous
statistical treatment of hardness ratios that properly deals with detected
photons as independent Poisson random variables and correctly deals with the
non-Gaussian nature of the error propagation. The method is Bayesian in nature,
and thus can be generalized to carry out a multitude of
source-population--based analyses. We verify our method with simulation
studies, and compare it with the classical method. We apply this method to real
world examples, such as the identification of candidate quiescent Low-mass
X-ray binaries in globular clusters, and tracking the time evolution of a flare
on a low-mass star.
