[ABRIDGED] The knowledge of the regular component of the Galactic magnetic field gives important information about the structure and dynamics of the Milky Way, as well as constitutes a basic tool to determine cosmic rays trajectories. It can also provide clear windows where primordial magnetic fields could be detected. We want 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 5-year WMAP data at 22 GHz. We have done a systematic study of a number of Galactic magnetic field models: axisymmetric, bisymmetric, 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 are computed. The model selection is done using a Bayesian approach. For each model, the posterior distributions are obtained and marginalised over the unwanted parameters to obtain the marginal 1-D probability distribution functions. In general, axisymmetric models provide a better description of the halo component, although attending to their goodness-of-fit, the rest of the models cannot be rejected. In the case of 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. 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 degrees, and a strong vertical field of 1 microG at z ~ 1 kpc. When a radial variation is fitted, models require fast variations.