Randomized benchmarking (RB) is the most commonly employed protocol for the characterization of unitary operations in quantum circuits due to its reasonable experimental requirements and robustness against state preparation and measurement (SPAM) errors. So far, the protocol has been limited to discrete or fermionic systems, whereas extensions to bosonic systems have been unclear for a long time due to challenges arising from the underlying bosonic Hilbert space. In this work, we close the gap for bosonic systems and develop an RB protocol to benchmark passive Gaussian transformations on any particle-number subspace, which we call passive bosonic RB. The protocol is built on top of the recently developed filtered RB framework [J. Helsen et al., PRX Quantum 3, 020357 (2022), M. Heinrich et al., Randomized benchmarking with random quantum circuits, arxiv:2212.06181 [quant-ph]] and is designed to isolate the multitude of exponential decays arising for passive bosonic transformations. We give explicit formulas and a Julia implementation for the necessary postprocessing of the experimental data. We also analyze the sampling complexity of passive bosonic RB by deriving analytical expressions for the variance. They show a mild scaling with the number of modes, suggesting that passive bosonic RB is experimentally feasible for a moderate number of modes. We focus on experimental settings involving Fock states and particle-number-resolving measurements, but also discuss Gaussian settings, deriving the first results for heterodyne measurements.