We present posterior distributions of parameters that characterize the nonperturbative initial input for the Balitsky-Kovchegov evolution equation. The BK equation evolves an initial dipole-target scattering amplitude at moderate $x_\mathrm{Bj}$ values toward smaller values. We use Bayesian inference to constrain the model parameters against the precise combined reduced cross section data from HERA, and obtain probability distributions for the free parameters. Our results show the data's preference for the anomalous dimension to be $\gamma \sim 1$. The distributions provide a rigorous method for propagating uncertainties of the BK initial condition for other CGC calculations. We demonstrate this for inclusive quark production in proton-proton collisions and nuclear modification factor for the total deep inelastic scattering cross section to be measured at the EIC.