We have revisited the problem of metallicity prediction of RR Lyrae stars from their near-infrared light curves in the Cousins I waveband. Our study is based on high-quality time-series photometry and state-of-the-art high-resolution spectroscopic abundance measurements of 80 fundamental-mode (RRab) and 24 first-overtone (RRc) stars, spanning ∼[−2.7,+0.18] dex and ∼[−3,−0.5] dex ranges, respectively. Employing machine-learning methods, we investigated various light-curve representations and regression models to identify their optimal form for our objective. Accurate new empirical relations between the [Fe/H] iron abundance and the light-curve parameters have been obtained using Bayesian regression for both RRab and RRc stars with mean absolute prediction errors of 0.16 and 0.18 dex, respectively. We found that earlier I-band [Fe/H] estimates had a systematic positive bias of up to ∼0.4 dex in the metal-poor regime. Our new predictive models were deployed on large ensembles of RR Lyrae stars to obtain photometric metallicity distribution functions (MDFs) for various old stellar populations in and around the Milky Way. We find that the mode of the old bulge component's MDF is approximately −1.4 dex, in remarkable agreement with the latest spectroscopic result. Furthermore, we derive MDF modes of −1.83, −2.13, and −1.77 dex for the Large and Small Magellanic Clouds, and the Sagittarius dwarf galaxy, respectively.