We propose a nonparametric approach for the identification of Wiener systems. We model the impulse response of the linear block and the static nonlinearity using Gaussian processes. The hyperparameters of the Gaussian processes are estimated using an iterative algorithm based on stochastic approximation expectation-maximization. In the iterations, we use elliptical slice sampling to approximate the posterior distribution of the impulse response and update the hyperparameter estimates. The same sampling is finally used to sample the posterior distribution and to compute point estimates. We compare the proposed approach with a parametric approach and a semi-parametric approach. In particular, we show that the proposed method has an advantage when a parametric model for the system is not readily available. (C) 2019 Elsevier Ltd. All rights reserved.