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.