  We propose a novel nonparametric regression framework subject to the positive
definiteness constraint. It offers a highly modular approach for estimating
covariance functions of stationary processes. Our method can impose positive
definiteness, as well as isotropy and monotonicity, on the estimators, and its
hyperparameters can be decided using cross validation. We define our estimators
by taking integral transforms of kernel-based distribution surrogates. We then
use the iterated density estimation evolutionary algorithm, a variant of
estimation of distribution algorithms, to fit the estimators. We also extend
our method to estimate covariance functions for point-referenced data. Compared
to alternative approaches, our method provides more reliable estimates for
long-range dependence. Several numerical studies are performed to demonstrate
the efficacy and performance of our method. Also, we illustrate our method
using precipitation data from the Spatial Interpolation Comparison 97 project.
