A stochastic eddy-damped quasi-normal Markovian (EDQNM) approach is used to investigate self-similar decaying isotropic turbulence at a high Reynolds number (400 ≤ Reλ ≤ 104). The realistic energy spectrum functional form recently proposed by Meyers & Menevau (Phys. Fluids, vol. 20, 2008, p. 065109) is generalized by considering some of the model constants as random parameters, since they escape measure in most experimental set-ups. The induced uncertainty on the solution is investigated, building response surfaces for decay power-law exponents of usual physical quantities. Large-scale uncertainties are considered, the emphasis being put on Saffman and Batchelor turbulences. The sensitivity of the solution to initial spectrum uncertainties is quantified through probability density functions of the decay exponents. It is observed that the initial spectrum shape at very large scales governs the long-time evolution, even at a high Reynolds number, a parameter which is not explicitly taken into account in many theoretical works. Therefore, a universal asymptotic behaviour in which kinetic energy decays as t−1 is not detected. However, this decay law is observed at finite Reynolds numbers with low probability for some initial conditions.