Summary: Upcoming galaxy redshift surveys promise to significantly improve current limits on primordial non-Gaussianity (PNG) through measurements of 2- and 3-point correlation functions in Fourier space. However, realizing the full potential of this dataset is contingent upon having both accurate theoretical models and optimized analysis methods. Focusing on the local model of PNG, parameterized by \(f_{\mathrm{NL}}\), we perform a Monte-Carlo Markov Chain analysis to confront perturbation theory predictions of the halo power spectrum and bispectrum in real space against a suite of N-body simulations. We model the halo bispectrum at tree-level, including all contributions linear and quadratic in \(f_{\mathrm{NL}}\), and the halo power spectrum at 1-loop, including tree-level terms up to quadratic order in \(f_{\mathrm{NL}}\), and all loops induced by local PNG linear in \(f_{\mathrm{NL}}\). Keeping the cosmological parameters fixed, we examine the effect of informative priors on the linear non-Gaussian bias parameter on the statistical inference of \(f_{\mathrm{NL}}\). A conservative analysis of the combined power spectrum and bispectrum, in which only loose priors are imposed and all parameters are marginalized over, can improve the constraint on by more than a factor of 5 relative to the power spectrum-only measurement. Imposing a strong prior on \(b_\phi\), or assuming bias relations for both \(b_\phi\) and \(b_{\phi\delta}\) (motivated by a universal mass function assumption), improves the constraints further by a factor of few. In this case, however, we find a significant systematic shift in the inferred value of \(f_{\mathrm{NL}}\) if the same range of wavenumber is used. Likewise, a Poisson noise assumption can lead to significant systematics, and it is thus essential to leave all the stochastic amplitudes free.