  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.
