The Milne–Eddington (M-E) atmosphere model is commonly adopted in the inversion of the magnetic fields in the solar photosphere. By applying the Levenberg–Marquardt algorithm or training a neural network (NN) model, the magnetic field vector can be quickly inferred from the Stokes profile but lacks reliable and statistically well-defined confidence intervals for parameters. To address this, we present an efficient Bayesian inference method called NNHMC, combining the NN model with the Hamiltonian Monte Carlo (HMC) algorithm. The NN model is used to speedily synthesize batches of synthetic Stokes profiles, accelerating the inference process. The HMC algorithm significantly improves sampling efficiency in high-dimensional parameter spaces and can handle large-scale data sets in batches. The spectropolarimetric observation of an active region obtained by the Hinode/spectropolarimeter (SP) is used to demonstrate the capability of the NNHMC method. The strength, inclination, and azimuth of the magnetic field and the line-of-sight velocity inferred with the NNHMC method are very similar to those derived with the MERLIN code. Furthermore, this study provided posterior distributions and uncertainties for these parameters. A test on the same hardware and software platform shows a speed increase of up to 2.5 orders of magnitude with respect to the traditional Markov Chain Monte Carlo method (without the NN, using the M-E atmosphere model), establishing the NNHMC method as a highly effective tool for Stokes inversion based on Bayesian inference.