A stochastic algorithm based on Metropolis Monte Carlo (MC) is presented for the size-extensive vibrational self-consistent field methods (XVSCF(n) and XVSCF[n]) for anharmonic molecular vibrations. The new MC-XVSCF methods substitute stochastic evaluations of a small number of high-dimensional integrals of functions of the potential energy surface (PES), which is sampled on demand, for diagrammatic equations involving high-order anharmonic force constants. This algorithm obviates the need to evaluate and store any high-dimensional partial derivatives of the potential and can be applied to the fully anharmonic PES without any Taylor-series approximation in an intrinsically parallelizable algorithm. The MC-XVSCF methods reproduce deterministic XVSCF calculations on the same Taylor-series PES in all energies, frequencies, and geometries. Calculations using the fully anharmonic PES evaluated on the fly with electronic structure methods report anharmonic effects on frequencies and geometries of much greater magnitude than deterministic XVSCF calculations, reflecting an underestimation of anharmonic effects in a Taylor-series approximation to the PES.