High-order sequential simulation techniques for complex and non-Gaussian spatially distributed variables were developed over the last few years. This paper presents a new extension where high-order statistics are inferred from the available hard data and then missing high-orders are borrowed from a training image. The inferred high-order statistics are then used to estimate spline coefficients that are subsequently employed to approximate conditional probability density functions as needed for the simulation process. The advantage of using orthogonal splines with respect to standard approaches is their ability to approximate a probability density function locally using not only high-order spatial cumulants for the whole range of values in data, but also partial cumulants calculated for particular ranges of values, such as extreme values. The proposed technique provides a general framework for simulation, both of continuous and categorical variables. Developments are tested on a synthetic data set.