Stokes inversion techniques, based on the analysis of the polarization in spectral lines, are the most powerful tools to obtain the information about the magnetic and thermodynamic quantities in the solar atmosphere. In this paper, we present a Milne-Eddington inversion code based on Bayesian inference, which is implemented with the Markov chain Monte Carlo simulation. This kind of simulation is a random sampling method to visit a point in the parameter space with a probability proportional to the distribution function, and furthermore provides a powerful way to infer the parameters of a model, their errors, and mutual correlations between each of them. Detailed tests of the code with synthetic profiles and Hinode/SP observations manifest the applicability of this code to infer physical quantities from Stokes profiles. Applying the code, we investigate the spatial distributions of physical quantities and their errors, and find that the errors of B, theta, phi, B-x, and B-y are smaller in the active regions than in the quiet regions. In contrast, the errors of B-z are much smaller in the quiet regions than in the plage and umbral regions. Besides, the study of the mutual correlations between each model parameter indicates that the thermodynamic parameters are strongly correlated, while the correlations between magnetic field parameters are relatively weak.