  We apply a suite of different estimators to the Quijote-PNG halo catalogues
to find the best approach to constrain Primordial non-Gaussianity (PNG) at
non-linear cosmological scales, up to $k_{\rm max} = 0.5 \, h\,{\rm Mpc}^{-1}$.
The set of summary statistics considered in our analysis includes the power
spectrum, bispectrum, halo mass function, marked power spectrum, and marked
modal bispectrum. Marked statistics are used here for the first time in the
context of PNG study. We perform a Fisher analysis to estimate their
cosmological information content, showing substantial improvements when marked
observables are added to the analysis. Starting from these summaries, we train
deep neural networks (NN) to perform likelihood-free inference of cosmological
and PNG parameters. We assess the performance of different subsets of summary
statistics; in the case of $f_\mathrm{NL}^\mathrm{equil}$, we find that a
combination of the power spectrum and a suitable marked power spectrum
outperforms the combination of power spectrum and bispectrum, the baseline
statistics usually employed in PNG analysis. A minimal pipeline to analyse the
statistics we identified can be implemented either with our ML algorithm or via
more traditional estimators, if these are deemed more reliable.
