  In this article, we introduce the frozen Gaussian sampling (FGS) algorithm to
solve the scalar wave equation in the high-frequency regime. The FGS algorithm
is a Monte Carlo sampling strategy based on the frozen Gaussian approximation,
which greatly reduces the computation workload in the wave propagation and
reconstruction. In this work, we propose feasible and detailed procedures to
implement the FGS algorithm to approximate scalar wave equations with Gaussian
initial conditions and WKB initial conditions respectively. For both initial
data cases, we rigorously analyze the error of applying this algorithm to wave
equations of dimensionality $d \geq 3$. In Gaussian initial data cases, we
prove that the sampling error due to the Monte Carlo method is independent of
the typical wave number. We also derive a quantitative bound of the sampling
error in WKB initial data cases. Finally, we validate the performance of the
FGS and the theoretical estimates about the sampling error through various
numerical examples, which include using the FGS to solve wave equations with
both Gaussian and WKB initial data of dimensionality $d = 1, 2$, and $3$.
