This paper presents an unsupervised algorithm for nonlinear unmixing of hyperspectral images. The proposed model assumes that the pixel reflectances result from a nonlinear function of the abundance vectors associated with the pure spectral components. We assume that the spectral signatures of the pure components and the nonlinear function are unknown. The first step of the proposed method estimates the abundance vectors for all the image pixels using a Bayesian approach an a Gaussian process latent variable model for the nonlinear function (relating the abundance vectors to the observations). The endmembers are subsequently estimated using Gaussian process regression. The performance of the unmixing strategy is first evaluated on synthetic data. The proposed method provides accurate abundance and endmember estimations when compared to other linear and nonlinear unmixing strategies. An interesting property is its robustness to the absence of pure pixels in the image. The analysis of a real hyperspectral image shows results that are in good agreement with state of the art unmixing strategies and with a recent classification method.