We present a comprehensive statistical analysis of the properties of Type Ia SN light curves in the near infrared using recent data from PAIRITEL and the literature. We construct a hierarchical Bayesian framework, incorporating several uncertainties
including photometric error, peculiar velocities, dust extinction and intrinsic variations, for coherent statistical inference. SN Ia light curve inferences are drawn from the global posterior probability of parameters describing both individual supernovae and the population conditioned on the entire SN Ia NIR dataset. The logical structure of the hierarchical model is represented by a directed acyclic graph. Fully Bayesian analysis of the model and data is enabled by an efficient MCMC algorithm exploiting the conditional structure using Gibbs sampling. We apply this framework to the JHK_s SN Ia light curve data. A new light curve model captures the observed J-band light curve shape variations. The intrinsic variances in peak absolute magnitudes are: sigma(M_J) = 0.17 +/- 0.03, sigma(M_H) = 0.11 +/- 0.03, and sigma(M_Ks) = 0.19 +/- 0.04. We describe the first quantitative evidence for correlations between the NIR absolute magnitudes and J-band light curve shapes, and demonstrate their utility for distance estimation. The average residual in the Hubble diagram for the training set SN at cz > 2000 km/s is 0.10 mag. The new application of bootstrap cross-validation to SN Ia light curve inference tests the sensitivity of the model fit to the finite sample and estimates the prediction error at 0.15 mag. These results demonstrate that SN Ia NIR light curves are as effective as optical light curves, and, because they are less vulnerable to dust absorption, they have great potential as precise and accurate cosmological distance indicators.
