  We apply Bayesian inference to determine the posterior likelihood
distribution for the parameters describing the initial condition of the
small-$x$ Balitsky-Kovchegov evolution equation at leading logarithmic
accuracy. The HERA structure function data is found to constrain most of the
model parameters well. In particular, we find that the HERA data prefers an
anomalous dimension $\gamma\approx 1$ unlike in previous fits where $\gamma>1$
which led to e.g. the unintegrated gluon distribution and the quark-target
cross sections not being positive definite. The determined posterior
distribution can be used to propagate the uncertainties in the non-perturbative
initial condition when calculating any other observable in the Color Glass
Condensate framework. We demonstrate this explicitly for the inclusive quark
production cross section in proton-proton collisions and by calculating
predictions for the nuclear modification factor for the $F_2$ structure
function in the EIC and LHeC/FCC-he kinematics.
