We perform population synthesis studies of different types of neutron stars (NSs) (thermally emitting isolated NSs, normal radio pulsars, magnetars) taking into account the magnetic field decay and using results from the most recent advances in NS cooling theory. For the first time, we confront our results with observations using simultaneously the log N–log S distribution for nearby isolated NSs, the log N–log L distribution for magnetars, and the distribution of radio pulsars in the diagram. For this purpose, we fix a baseline NS model (all microphysics input), and other relevant parameters to standard values (velocity distribution, mass spectrum, birth rates, etc.), allowing us to vary the initial magnetic field strength. We find that our theoretical model is consistent with all sets of data if the initial magnetic field distribution function follows a lognormal law with 〈log (B0/G)〉∼ 13.25 and . The typical scenario includes about 10 per cent of NSs born as magnetars, significant magnetic field decay during the first million years of a NS life (only about a factor of 2 for low-field NSs but more than an order of magnitude for magnetars), and a mass distribution function dominated by low-mass objects. This model explains satisfactorily all known populations. Evolutionary links between different subclasses may exist, although robust conclusions are not yet possible.