Different viscoelastic models and characteristics are commonly used to describe, analyze, compare and improve the mechanical properties of polymers. A time-dependent linear relaxation modulus next to frequency-domain storage and loss moduli are the basic rheological material functions of polymers. The exponential Maxwell model and the exponential stretched Kohlrausch–Williams–Watts model are, probably, the most known linear rheological models of polymers. There are different identification methods for such models, some of which are dedicated to specific models, while others are general in nature. However, the identification result, i.e., the best model, always depends on the specific experimental data on the basis of which it was determined. When the rheological stress relaxation test is performed, the data are composed of the sampling instants used in the test and on the measurements of the relaxation modulus of the real material. To build a relaxation modulus model that does not depend on sampling instants is a fundamental concern. The problem of weighted least-squares approximation of the real relaxation modulus is discussed when only the noise-corrupted time-measurements of the relaxation modulus are accessible for identification. A wide class of models, that are continuous, differentiable and Lipschitz with respect to parameters, is considered for the relaxation modulus approximation. The main results concern the models that are selected asymptotically as the number of measurements tends to infinity. It is shown that even when the true relaxation modulus description is completely unknown, the approximate optimal model parameters can be derived from the measurement data that are obtained for sampling instants that are selected randomly due to the appropriate randomization introduced whenever certain conditions regarding the adopted class of models are satisfied. It is shown that the most commonly used stress relaxation models, the Maxwell and Kohlrausch–Williams–Watts models, satisfy these conditions. Since the practical problems of the identification of relaxation modulus models are usually ill posed, Tikhonov regularization is applied to guarantee the stability of the regularized solutions. The approximate optimal model is a strongly consistent estimate of the regularized model that is optimal in the sense of the deterministic integral weighted square error. An identification algorithm leading to the best regularized model is presented. The stochastic-type convergence analysis is conducted for noise-corrupted relaxation modulus measurements, and the exponential convergence rate is proved. Numerical studies for different models of the relaxation modulus used in the polymer rheology are presented for the material described by a bimodal Gauss-like relaxation spectrum. Numerical studies have shown that if appropriate randomization is introduced in the selection of sampling instants, then optimal regularized models of the relaxation modulus being asymptotically independent of these time instants can be recovered from the stress relaxation experiment data. The robustness of the identification algorithm to measurement noises was demonstrated both by analytical and numerical analyses.