Toward a robust inference method for the galaxy bispectrum: likelihood function and model selection