Context. The knowledge of the regular (large scale) component of the Galactic magnetic field gives important information about the structure and dynamics of the Milky Way, and also constitutes a basic tool to determine cosmic ray trajectories. It can also provide clear windows where primordial magnetic fields could be detected. Aims: We aim to obtain the regular (large scale) pattern of the magnetic field distribution of the Milky Way that better fits the polarized synchrotron emission as seen by the WMAP satellite in the 5 years data at 22 GHz. Methods: We have done a systematic study of a number of Galactic magnetic field models: axisymmetric (with and without radial dependence on the field strength), bisymmetric (with and without radial dependence), logarithmic spiral arms, concentric circular rings with reversals and bi-toroidal. We have explored the parameter space defining each of these models using a grid-based approach. In total, more than one million models were computed. The model selection was done using a Bayesian approach. For each model, the posterior distributions were obtained and marginalized over the unwanted parameters to obtain the marginal (one-parameter) probability distribution functions. Results: In general, axisymmetric models provide a better description of the halo component, although with regard to their goodness-of-fit, the other models cannot be rejected. In the case of the disk component, the analysis is not very sensitive for obtaining the disk large-scale structure, because of the effective available area (less than 8% of the whole map and less than 40% of the disk). Nevertheless, within a given family of models, the best-fit parameters are compatible with those found in the literature. Conclusions: The family of models that better describes the polarized synchrotron halo emission is the axisymmetric one, with magnetic spiral arms with a pitch angle of ≈24°, and a strong vertical field of 1 μG at z ≈ 1 kpc. When a radial variation is fitted, models require fast variations.