  In the new precision era for cosmic ray astrophysics, theoretical predictions
cannot content themselves with average trends, but need to correctly take into
account intrinsic uncertainties. The space-time discreteness of the cosmic ray
sources, together with a substantial ignorance of their precise epochs and
locations (with the possible exception of the most recent and close ones) play
an important role in this sense. We elaborate a statistical theory to deal with
this problem, relating the composite probability P({\Psi}) to obtain a flux
{\Psi} at the Earth and the single-source probability p({\psi}) to contribute
with a flux {\psi}. The main difficulty arises from the fact that p({\psi}) is
a fat tail distribution, characterized by power-law or broken power-law
behavior up to very large fluxes for which central limit theorem does not hold,
and leading to well-known stable laws as opposed Gaussian. We find that
relatively simple recipes provide a satisfactory description of the probability
P({\Psi}). We also find that a naive Gaussian fit to simulation results would
underestimate the probability of very large fluxes, that is, several times
above the average, while overestimating the probability of relatively milder
excursions. At large energies, large flux fluctuations are prevented by causal
considerations, while at low energies a partial knowledge on the recent and
nearby population of sources plays an important role. A few proposal have been
discussed in the literature to account for spectral breaks recently reported in
cosmic ray data in terms of local contributions. We apply our newly developed
theory to assess their probabilities, finding that they are relatively small.
