We present a new method based on a Bayesian hierarchical model to extract constraints on cosmological parameters from SNIa data obtained with the SALT-II lightcurve fitter. We demonstrate with simulated data sets that our method delivers tighter statistical constraints on the cosmological parameters over 90% of the time, that it reduces statistical bias typically by a factor ~ 2-3 and that it has better coverage properties than the usual chi-squared approach. As a further benefit, a full posterior probability distribution for the dispersion of the intrinsic magnitude of SNe is obtained. We apply this method to recent SNIa data, and by combining them with CMB and BAO data we obtain Omega_m=0.28 +/- 0.02, Omega_Lambda=0.73 +/- 0.01 (assuming w=-1) and Omega_m=0.28 +/- 0.01, w=-0.90 +/- 0.05 (assuming flatness; statistical uncertainties only). We constrain the intrinsic dispersion of the B-band magnitude of the SNIa population, obtaining sigma_mu^int = 0.13 +/- 0.01 [mag]. Applications to systematic uncertainties will be discussed in a forthcoming paper.