We propose a new stochastic algorithm for computing useful Bayesian estimators of hidden Markov random field (HMRF) models that we call exploration/selection/estimation (ESE) procedure. The algorithm is based on an optimization algorithm of O. François, called the exploration/selection (E/S) algorithm. The novelty consists of using the a posteriori distribution of the HMRF, as exploration distribution in the E/S algorithm. The ESE procedure computes the estimation of the likelihood parameters and the optimal number of region classes, according to global constraints, as well as the segmentation of the image. In our formulation, the total number of region classes is fixed, but classes are allowed or disallowed dynamically. This framework replaces the mechanism of the split-and-merge of regions that can be used in the context of image segmentation. The procedure is applied to the estimation of a HMRF color model for images, whose likelihood is based on multivariate distributions, with each component following a Beta distribution. Meanwhile, a method for computing the maximum likelihood estimators of Beta distributions is presented. Experimental results performed on 100 natural images are reported. We also include a proof of convergence of the E/S algorithm in the case of nonsymmetric exploration graphs.