A simple and numerically efficient approach to Wigner transforms and classical Liouville dynamics in phase space is presented. (1) The Wigner transform can be obtained with a given accuracy by optimal decomposition of an initial quantum-mechanical wave function in terms of a minimal set of Gaussian wave packets. (2) The solution of the classical Liouville equation within the locally quadratic approximation of the potential energy function requires a representation of the density in terms of an ensemble of narrow Gaussian phase-space packets. The corresponding equations of motion can be efficiently solved by a modified leap-frog integrator. For both problems the use of Monte Carlo based techniques allows numerical calculation in multidimensional cases where grid-based methods such as fast Fourier transforms are prohibitive. In total, the proposed strategy provides a practical and efficient tool for classical Liouville dynamics with quantum-mechanical initial states.