Stokes inversion techniques are very powerful methods for obtaining information on the thermodynamic and magnetic properties of solar and stellar atmospheres. Most of the existing inversion codes are designed for finding the optimum solution to the nonlinear inverse problem. However, to obtain the location of potentially multimodal solutions, degeneracies, and the uncertainties of each parameter from the inversions, algorithms such as Markov chain Monte Carlo require to evaluate the model thousand of times. Variational methods are a quick alternative by approximating the posterior distribution by a parametrized distribution. In this study, we explore a highly flexible variational method, known as normalizing flows, to return Bayesian posterior probabilities for solar observations. We illustrate the ability of the method using a simple Milne-Eddington model and a complex non-LTE inversion. The training procedure need only be performed once for a given prior parameter space and the resulting network can then generate samples describing the posterior distribution several orders of magnitude faster than existing techniques.