The halo occupation distribution (HOD) framework is an empirical method to describe the connection between dark matter halos and galaxies, which is constrained by small scale clustering data. Efficient fitting procedures are required to scan the HOD parameter space. This paper describes such a method based on Gaussian Processes to iteratively build a surrogate model of the posterior of the likelihood surface from a reasonable amount of likelihood computations, typically two orders of magnitude less than standard Monte Carlo Markov chain algorithms. Errors in the likelihood computation due to stochastic HOD modelling are also accounted for in the method we propose. We report results of reproducibility, accuracy and stability tests of the method derived from simulation, taking as a test case star-forming emission line galaxies, which constitute the main tracer of the Dark Energy Spectroscopic Instrument and have so far a poorly constrained galaxy-halo connection from observational data.