We present a method, based on Bayesian statistics, to fit the dust emission parameters in the far-infrared and submillimeter wavelengths. The method estimates the dust temperature and spectral emissivity index, plus their relationship, taking into account properly the statistical and systematic uncertainties. We test it on three sets of simulated sources detectable by the Herschel Space Observatory in the PACS and SPIRE spectral bands (70-500 micron), spanning over a wide range of dust temperatures. The simulated observations are a one-component Interstellar Medium, and two two-component sources, both warm (HII regions) and cold (cold clumps). We first define a procedure to identify the better model, then we recover the parameters of the model and measure their physical correlations by means of a Monte Carlo Markov Chain algorithm adopting multi-variate Gaussian priors. In this process we assess the reliability of the model recovery, and of parameters estimation. We conclude that the model and parameters are properly recovered only under certain circumstances, and that false models may be derived in some case. We applied the method to a set of 91 starless cold clumps in an inter-arm region of the Galactic Plane with low star formation activity, observed by Herschel in the Hi-GAL survey. Our results are consistent with a temperature independent spectral index.