This paper develops an algorithmic approach for obtaining estimates of the weight enumerators of Reed-Muller (RM) codes. Our algorithm is based on a technique for estimating the partition functions of spin systems, which in turn employs a sampler that produces codewords according to a suitably defined Gibbs distribution. We apply our method to moderate-blocklength RM codes and derive approximate values of their weight enumerators. We observe that the rates of the weight enumerator estimates returned by our method are close to the true rates when these rates are either known or computable by brute-force search; in other cases, our computations provide provably robust estimates. As a byproduct, our sampling algorithm also allows us to obtain estimates of the weight spectra of RM codes. We illustrate our methods by providing estimates of the hitherto unknown weight enumerators of the RM code and the exact weight spectra of the RM and RM codes.