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 statitical bias typically by a factor \sim 2 - 3 and that it has better coverage properties than the usual \chi ^ { 2 } 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 \pm 0.02 , \Omega _ { \Lambda } = 0.73 \pm 0.01 ( assuming w = -1 ) and \Omega _ { m } = 0.28 \pm 0.01 ,w = -0.90 \pm 0.05 ( assuming flatness ; statistical uncertainties only ) . We constrain the intrinsic dispersion of the B-band magnitude of the SNIa population , obtaining \sigma _ { \mu } ^ { \text { int } } = 0.13 \pm 0.01 [ mag ] . Applications to systematic uncertainties will be discussed in a forthcoming paper .