Considerable development has been observed in the area of applying fractional-order rheological models to describe the viscoelastic properties of miscellaneous materials in the last few decades together with the increasingly stronger adoption of fractional calculus. The fractional Maxwell model is the best-known non-integer-order rheological model. A weighted least-square approximation problem of the relaxation modulus by the fractional Maxwell model is considered when only the time measurements of the relaxation modulus corrupted by additive noises are accessible for identification. This study was dedicated to the determination of the model, optimal in the sense of the integral square weighted model quality index, which does not depend on the particular sampling points applied in the stress relaxation experiment. It is proved that even when the real description of the material relaxation modulus is entirely unknown, the optimal fractional Maxwell model parameters can be recovered from the relaxation modulus measurements recorded for sampling time points selected randomly according to respective randomization. The identified model is a strongly consistent estimate of the desired optimal model. The exponential convergence rate is demonstrated both by the stochastic convergence analysis and by the numerical studies. A simple scheme for the optimal model identification is given. Numerical studies are presented for the materials described by the short relaxation times of the unimodal Gauss-like relaxation spectrum and the long relaxation times of the Baumgaertel, Schausberger and Winter spectrum. These studies have shown that the appropriate randomization introduced in the selection of sampling points guarantees that the sequence of the optimal fractional Maxwell model parameters asymptotically converge to parameters independent of these sampling points. The robustness of the identified model to the measurement disturbances was demonstrated by analytical analysis and numerical studies.