Preferential sampling of high grade zones is a common practice in the drilling campaigns of ore deposits. This may lead to a bias in the statistics derived from these clustered datasets. The global mean and variance, particularly, tend to be overestimated and high-order statistics, such as the skewness and kurtosis are distorted. High-order simulation uses high-order spatial statistics, known as cumulants, to approximate non-Gaussian local distributions of grades conditioned by neighbouring samples. Data clustering affects the inference of the spatial cumulants, and consequently, the conditional distributions fitted using them. A weighted cumulant estimator is proposed to account for data clustering. This estimator is implemented in the program HOSC, used for calculating the experimental spatial cumulants of irregularly spaced datasets. Additionally, the fitting of non-Gaussian conditional distributions using Legendre polynomial series derived from weighted spatial cumulants is proposed and implemented in the algorithm for simulation with high-order statistics, HOSIM. The geological model and diamond drilling dataset corresponding to the Apensu gold deposit in Ghana is used to illustrate these implementations and to compare the new results with the realizations and statistics produced without accounting for data clustering. This example shows that as expected, the data cdf and its statistics are reproduced by simulation with declustered high-order statistics.