Primordial non-Gaussianity introduces a scale-dependent variation in the clustering of density peaks corresponding to rare objects. This variation, parametrized by the bias, is investigated on scales where a linear perturbation theory is sufficiently accurate. The bias is obtained directly in real space by comparing the one-and two-point probability distributions of density fluctuations. We show that these distributions can be reconstructed using a bivariate Edgeworth series, presented here up to an arbitrarily high order. The Edgeworth formalism is shown to be well-suited for 'local' cubic-order non-Gaussianity parametrized by g NL . We show that a strong scale-dependence in the bias can be produced by g NL of order 10 5 , consistent with CMB constraints. On correlation length of ∼ 100 Mpc, current constraints on g NL still allow the bias for the most massive clusters to be enhanced by 20 − 30% of the Gaussian value. We further examine the bias as a function of mass scale, and also explore the relationship between the clustering and the abundance of massive clusters in the presence of g NL . We explain why the Edgeworth formalism, though technically challenging, is a very powerful technique for constraining high-order non-Gaussianity with large-scale structures. * Electronic address: siri@astro.ox.ac.uk 1