  A method based on Monte Carlo techniques is presented for evaluating
thermonuclear reaction rates. We begin by reviewing commonly applied procedures
and point out that reaction rates that have been reported up to now in the
literature have no rigorous statistical meaning. Subsequently, we associate
each nuclear physics quantity entering in the calculation of reaction rates
with a specific probability density function, including Gaussian, lognormal and
chi-squared distributions. Based on these probability density functions the
total reaction rate is randomly sampled many times until the required
statistical precision is achieved. This procedure results in a median (Monte
Carlo) rate which agrees under certain conditions with the commonly reported
recommended "classical" rate. In addition, we present at each temperature a low
rate and a high rate, corresponding to the 0.16 and 0.84 quantiles of the
cumulative reaction rate distribution. These quantities are in general
different from the statistically meaningless "minimum" (or "lower limit") and
"maximum" (or "upper limit") reaction rates which are commonly reported.
Furthermore, we approximate the output reaction rate probability density
function by a lognormal distribution and present, at each temperature, the
lognormal parameters miu and sigma. The values of these quantities will be
crucial for future Monte Carlo nucleosynthesis studies. Our new reaction rates,
appropriate for bare nuclei in the laboratory, are tabulated in the second
paper of this series (Paper II). The nuclear physics input used to derive our
reaction rates is presented in the third paper of this series (Paper III). In
the fourth paper of this series (Paper IV) we compare our new reaction rates to
previous results.
