  In this work, we study the semi-classical limit of the Schr\"odinger equation
with random inputs, and show that the semi-classical Schr\"odinger equation
produces $O(\varepsilon)$ oscillations in the random variable space. With the
Gaussian wave packet transform, the original Schr\"odinger 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 $\varepsilon$ dependent oscillations, 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,
i.e. simulating the $w$ equation, it is sufficient to use $\varepsilon$
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.
