We present a detailed description of the recently proposed numerical renormalization group method for models of quantum impurities coupled to a bosonic bath. Specifically, the method is applied to the spin-boson model, both in the Ohmic and sub-Ohmic cases. We present various results for static as well as dynamic quantities and discuss details of the numerical implementation, e.g., the discretization of a bosonic bath with arbitrary continuous spectral density, the suitable choice of a finite basis in the bosonic Hilbert space, and questions of convergence with respect to truncation parameters. The method is shown to provide high-accuracy data over the whole range of model parameters and temperatures, which are in agreement with exact results and other numerical data from the literature.