The classical Wigner model is one way to approximate the quantum dynamics of atomic nuclei. Here, a new method is presented for sampling the initial quantum mechanical distribution that is required in the classical Wigner model. The new method is tested for the position, position-squared, momentum, and momentum-squared autocorrelation functions for a one-dimensional quartic oscillator and double well potential as well as a quartic oscillator coupled to harmonic baths of different sizes. Two versions of the new method are tested and shown to possibly be useful. Both versions always converge toward the classical Wigner limit. For the one-dimensional cases, some results that are essentially converged to the classical Wigner limit are acquired and others are not far off. For the multi-dimensional systems, the convergence is slower, but approximating the sampling of the harmonic bath with classical mechanics was found to greatly improve the numerical performance. For the double well, the new method is noticeably better than the Feynman-Kleinert linearized path integral method at reproducing the exact classical Wigner results, but they are equally good at reproducing exact quantum mechanics. The new method is suggested as being interesting for future tests on other correlation functions and systems.