We investigated the structure of the low density regions of the inner crust of neutron stars using the Hartree-Fock-Bogoliubov (HFB) model to predict the proton content Z of the nuclear clusters and, together with the lattice spacing, the proton content of the crust as a function of the total baryonic density ρ b . The exploration of the energy surface in the (Z, ρ b ) configuration space and the search for the local minima require thousands of calculations. Each of them implies an HFB calculation in a box with a large number of particles, thus making the whole process very demanding. In this work, we apply a statistical model based on a Gaussian Process Emulator that makes the exploration of the energy surface ten times faster. We also present a novel treatment of the HFB equations that leads to an uncertainty on the total energy of ≈ 4 keV per particle. Such a high precision is necessary to distinguish neighbour configurations around the energy minima.