In this paper we propose a method, based on a Set Membership approach, for the estimation of nonlinear regressions models. At the contrary of most of the existing identification approaches, the method presented in this paper does not need any assumption about the functional form of the model to be identified, but uses only some prior information on its regularity and on the size of noise corrupting the measurements. The aim is to evaluate not only a nominal model but a model set, describing the inherent uncertainty of the regression function coming from finite and noise corrupted data. This is obtained by computing the optimal bounds on the regression function , i.e. its tightest lower and upper bounds compatible with measured data and with the given assumptions on the regression function and on noise. Moreover, necessary and sufficient conditions are given for validating the prior assumptions. The effectiveness of the method is tested on a water heater identification problem, where the obtained models are compared in simulation with other nonlinear models obtained by neural networks, Just In Time and Fuzzy approaches.