  We analyze a set of 89 Type Ia supernovae (SN Ia) that have both optical and
near-infrared (NIR) photometry to derive distances and construct low redshift
($z < 0.04$) Hubble diagrams. We construct mean light curve (LC) templates
using a hierarchical Bayesian model. We explore both Gaussian process (GP) and
template methods for fitting the LCs and estimating distances, while including
peculiar velocity and photometric uncertainties. For the 56 SN Ia with both
optical and NIR observations near maximum light, the GP method yields a
NIR-only Hubble-diagram with a RMS of $0.117 \pm 0.014$ mag when referenced to
the NIR maxima. For each NIR band, a comparable GP method RMS is obtained when
referencing to NIR-max or B-max. Using NIR LC templates referenced to B-max
yields a larger RMS value of $0.138 \pm 0.014$ mag. Fitting the corresponding
optical data using standard LC fitters that use LC shape and color corrections
yields larger RMS values of $0.179 \pm 0.018$ mag with SALT2 and $0.174 \pm
0.021$ mag with SNooPy. Applying our GP method to subsets of SN Ia NIR LCs at
NIR maximum light, even without corrections for LC shape, color, or host-galaxy
dust reddening, provides smaller RMS in the inferred distances, at the $\sim
2.3 - 4.1\sigma$ level, than standard optical methods that do correct for those
effects. Our ongoing RAISIN program on the Hubble Space Telescope will exploit
this promising infrared approach to limit systematic errors when measuring the
expansion history of the universe to constrain dark energy.
