We describe the generalization of Wilson’s numerical renormalization group method to quantum impurity models with a bosonic bath, providing a general nonperturbative approach to bosonic impurity models which can access exponentially small energies and temperatures. As an application, we consider the spin-boson model, describing a two-level system coupled to a bosonic bath with power-law spectral density, [Formula presented]. We find clear evidence for a line of continuous quantum phase transitions for sub-Ohmic bath exponents [Formula presented]; the line terminates in the well-known Kosterlitz-Thouless transition at [Formula presented]. Contact is made with results from perturbative renormalization group, and various other applications are outlined.