  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,000, 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.
