We model the evolution of active galactic nuclei by constructing their radio luminosity functions. We use a set of surveys of varying area and depth, namely the deep COSMOS survey of AGN sources, the wide shallow 3CRR, 7C and 6CE surveys, containing together AGNs, and the intermediate XXL-North and South fields consisting of and sources, respectively. We also used the CENSORS, BRL, Wall Peacock and Config surveys, consisting respectively of , , and sources. Together, these surveys numbered AGN sources and constrained the luminosity functions at high redshift and over a wide range of luminosities (up to and . We concentrate on parametric methods within the Bayesian framework and show that the luminosity-dependent density evolution (LDDE) model fits the data best, with evidence ratios varying from "strong" () to "decisive" () according to the Jeffreys interpretation. We determine the number density, luminosity density and kinetic luminosity density as a function of redshift, and observe a flattening of these functions at higher redshifts, not present in simpler models, which we explain by our use of the LDDE model. Finally, we divide our sample into subsets according to the stellar mass of the host galaxies in order to investigate a possible bimodality in evolution. We found a difference in LF shape and evolution between these subsets. All together, these findings point to a physical picture where the evolution and density of AGN cannot be explained well by simple models but require more complex models either via AGN sub-populations where the total AGN sample is divided into subsamples according to various properties such as, for example, optical properties and stellar mass, or via luminosity-dependent functions.