We show how importance sampling can be used to reconstruct the statistics of
rare cosmological fluctuations in stochastic inflation. We have developed a
publicly available package, PyFPT, that solves the first-passage time problem
of generic one-dimensional Langevin processes. In the stochastic-$\delta N$
formalism, these are related to the curvature perturbation at the end of
inflation. We apply this method to quadratic inflation, where the existence of
semi-analytical results allows us to benchmark our approach. We find excellent
agreement within the estimated statistical error, both in the drift- and
diffusion-dominated regimes. The computation takes at most a few hours on a
single CPU, and can reach probability values corresponding to less than one
Hubble patch per observable universe at the end of inflation. With direct
sampling, this would take more than the age of the universe to simulate even
with the best current supercomputers. As an application, we study how the
presence of large-field boundaries might affect the tail of the probability
distribution. We also find that non-perturbative deviations from Gaussianity
are not always of the simple exponential type.
