We propose a probabilistic approach for estimating parameters of an option pricing model from a set of observed option prices. Our approach is based on a stochastic optimization algorithm which generates a random sample from the set of global minima of the in-sample pricing error and allows for the existence of multiple global minima. Starting from an IID population of candidate solutions drawn from a prior distribution of the set of model parameters, the population of parameters is updated through cycles of independent random moves followed by "selection" according to pricing performance. We examine conditions under which such an evolving population converges to a sample of calibrated models. The heterogeneity of the obtained sample can then be used to quantify the degree of ill-posedness of the inverse problem: it provides a natural example of a coherent measure of risk, which is compatible with observed prices of benchmark ("vanilla") options and takes into account the model uncertainty resulting from incomplete identification of the model. We describe in detail the algorithm in the case of a diffusion model, where one aims at retrieving the unknown local volatility surface from a finite set of option prices, and illustrate its performance on simulated and empirical data sets of index options.