We use Bayesian inference to constrain the parameters describing the initial amplitude input to the Balitsky-Kovchegov evolution equation at next-to-leading order accuracy against precise HERA total inclusive cross section and heavy quark data. The datasets are found to provide stringent constraints and, with consistent NLO treatment, a successful description of the data is obtained. The posterior distributions define the theoretical uncertainites that surround the non-perturbative initial condition and, thus, provide a way to propagate said uncertainties to CGC calculations at NLO.