Precision long-baseline neutrino experiments use near-detector data to constrain systematic uncertainties on the unoscillated neutrino flux, a prerequisite for accurate oscillation-parameter measurements at the far detector. When the constrained likelihood is high-dimensional and non-Gaussian, this procedure demands advanced statistical treatment. Here we show that normalizing flows provide faithful and portable likelihood models for this problem. Leveraging an initial Gaussian approximation of the likelihood, we train a hybrid architecture combining coupling transformations and autoregressive spline flows. We demonstrate the method on a representative neardetector likelihood replica with 110 systematic uncertainty parameters, 10 of which explicitly introduce non-Gaussianities in the constrained likelihood. The trained model achieves a relative effective sample size of 98%, compared with about 5% for the Gaussian approximation, and reproduces a Markov Chain Monte Carlo reference while remaining closed-form, samplable, and pointwise evaluable, making it suitable to uncertainty propagation in downstream physics analyses.