Nonlinear energy sinks (NES) are a promising technique to achieve vibration mitigation. Through nonlinear stiffness properties, NES are able to passively and irreversibly absorb energy. Unlike the traditional Tuned Mass Damper (TMD), NES absorb energy from a wide range of frequencies. Many studies have focused on NES behavior and dynamics, but few have addressed the optimal design of NES. Design considerations of NES are of prime importance as it has been shown that NES dynamics exhibit an acute sensitivity to uncertainties. In fact, the sensitivity is so marked that NES efficiency is near-discontinuous and can switch from a high to a low value for a small perturbation in design parameters or loading conditions. This article presents an approach for the probabilistic design of NES which accounts for random design and aleatory variables as well as response discontinuities. In order to maximize the mean efficiency, the algorithm is based on the identification of regions of the design and aleatory space corresponding to markedly different NES efficiencies. This is done through a sequence of approximated sub-problems constructed from clustering, Kriging approximations, a support vector machine, and Monte-Carlo simulations. The refinement of the surrogates is performed locally using a generalized max-min sampling scheme which accounts for the distributions of random variables. The sampling scheme also makes use of the predicted variance of the Kriging surrogates for the selection of aleatory variables values. The proposed algorithm is applied to three example problems of varying dimensionality, all including an aleatory excitation applied to the main system. The stochastic optima are compared to NES optimized deterministically.