In this work, we study the semiclassical limit of the Schrödinger equation with random inputs, and show that the semiclassical Schrödinger equation produces O(epsilon) oscillations in the z variable in general. However, with the Gaussian wave packet transform, the original Schrödinger equation is mapped to an ODE system for the wave packet parameters coupled with a PDE for the quantity w in rescaled variables. Further, we show that the w equation does not produce epsilon dependent oscillations in the rescaled spatial variable, and thus it is more amenable for numerical simulations. We propose multi-level sampling strategy in implementing the Gaussian wave packet transform, where in the most costly part, simulating the w equation, it is sufficient to use ε independent samples. We also provide extensive numerical tests as well as meaningful numerical experiments to justify the properties of the numerical algorithm, and hopefully shed light on possible future directions.