Monte Carlo (MC) methods have been found useful in estimation of variance parameters for large data and complex models with many variance components (VC), with respect to both computer memory and computing time. A disadvantage has been a fluctuation in round-to-round values of estimates that makes the estimation of convergence challenging. Furthermore, with Newton-type algorithms, the approximate Hessian matrix might have sufficient accuracy, but the inaccuracy in the gradient vector exaggerates the round-to-round fluctuation to intolerable. In this study, the reuse of the same random numbers within each MC sample was used to remove the MC fluctuation. Simulated data with six VC parameters were analysed by four different MC REML methods: expectation-maximization (EM), Newton-Raphson (NR), average information (AI) and Broyden's method (BM). In addition, field data with 96 VC parameters were analysed by MC EM REML. In all the analyses with reused samples, the MC fluctuations disappeared, but the final estimates by the MC REML methods differed from the analytically calculated values more than expected especially when the number of MC samples was small. The difference depended on the random numbers generated, and based on repeated MC AI REML analyses, the VC estimates were on average non-biased. The advantage of reusing MC samples is more apparent in the NR-type algorithms. Smooth convergence opens the possibility to use the fast converging Newton-type algorithms. However, a disadvantage from reusing MC samples is a possible "bias" in the estimates. To attain acceptable accuracy, sufficient number of MC samples need to be generated.