Jets from supermassive black holes in the centers of active galaxies are the most powerful persistent sources of electromagnetic radiation in the Universe. Jets emit non-thermal radiation, with luminosity as high as 1048 erg/s, and their spectral energy distribution (SED) spans many decades in photon energy, i.e. from radio waves to gamma rays. To infer the physical conditions in the otherwise out-of-reach regions of extragalactic jets we usually rely on SED fitting, i.e. we compare observed SEDs with radiative models in order to measure the goodness-of-fit of the model to the data, and determine the best-fit values. Radiative models for jet emission are primarily numerical, as they require the solution of a system of stiff coupled partial differential equations describing the evolution of the distribution functions of radiating particles (relativistic leptons and/or hadrons) and photons. Such numerical codes tend to have a high computational complexity, especially in the case of hadronic models; each run can last from a few minutes up to a couple of hours. If such a model were to be used in a Markov Chain Monte Carlo (MCMC) algorithm, the overall execution time would become prohibitively long. In this work, machine learning is used to tackle the problem of high computational complexity in order to reduce the model evaluation time. An existing numerical code is executed in order to compose a leptonic dataset, which in turn is used to train a neural network. The trained neural network can make predictions in a few milliseconds and thus can be used instead of the actual numerical code. We use the trained network to fit observational data from an extra-galactic jet using Bayesian methods based on random walks and nested sampling algorithms. Our results demonstrate that our approach offers a viable alternative that could potentially replace expensive algorithms, especially when searching multi-parameter spaces.