This is the first paper in a series devoted to systematic study of the size and structure of the broad-line region (BLR) in active galactic nuclei (AGNs) using reverberation mapping (RM) data. We employ a recently developed Bayesian approach that statistically describes the variabibility as a damped random walk process and delineates the BLR structure using a flexible disk geometry that can account for a variety of shapes, including disks, rings, shells, and spheres. We allow for the possibility that the line emission may respond non-linearly to the continuum, and we detrend the light curves when there is clear evidence for secular variation. We use a Markov Chain Monte Carlo implementation based on Bayesian statistics to recover the parameters and uncertainties for the BLR model. The corresponding transfer function is obtained self-consistently. We tentatively constrain the virial factor used to estimate black hole masses; more accurate determinations will have to await velocity-resolved RM data. Application of our method to RM data with Hbeta monitoring for about 40 objects shows that the assumed BLR geometry can reproduce quite well the observed emission-line fluxes from the continuum light curves. We find that the Hbeta BLR sizes obtained from our method are on average ~20% larger than those derived from the traditional cross-correlation method. Nevertheless, we still find a tight BLR size-luminosity relation with a slope of alpha=0.55±0.03 and an intrinsic scatter of ~0.18 dex. In particular, we demonstrate that our approach yields appropriate BLR sizes for some objects (such as Mrk 142 and PG 2130+099) where traditional methods previously encountered difficulties. ; Comment: 17 pages, 10 figures, 2 tables; minor reversion to match the published version