The current and upcoming generations of gravitational wave experiments represent an exciting step forward in terms of detector sensitivity and performance. For example, key upgrades at the LIGO, Virgo and KAGRA facilities will see the next observing run (O4) probe a spatial volume around four times larger than the previous run (O3), and design implementations for, e.g., the Einstein Telescope, Cosmic Explorer, and LISA experiments are taking shape to explore a wider frequency range and probe cosmic distances. In this context, however, a number of very real data analysis problems face the gravitational wave community. For example, it will be critical to develop tools and strategies to analyze (among other scenarios) signals that arrive coincidentally in detectors, longer signals that are in the presence of nonstationary noise or other shorter transients, as well as noisy, potentially correlated, coherent stochastic backgrounds. With these challenges in mind, we develop peregrine, a new sequential simulation-based inference approach designed to study broad classes of gravitational wave signal. In this work, we describe the method and implementation, before demonstrating its accuracy and robustness through direct comparison with established likelihood-based methods. Specifically, we show that we are able to fully reconstruct the posterior distributions for every parameter of a spinning, precessing compact binary coalescence using one of the most physically detailed and computationally expensive waveform approximants (SEOBNRv4PHM). Crucially, we are able to do this using only 2% of the waveform evaluations that are required in, e.g., nested sampling approaches. Finally, we provide some outlook as to how this level of simulation efficiency and flexibility in the statistical analysis could allow peregrine to tackle these current and future gravitational wave data analysis problems