The Extragalactic Background Light (EBL) is the main radiation field responsible for
attenuating extragalactic gamma-ray emission at very high energies, but its precise spectral
intensity is not fully determined. Therefore, disentangling propagation effects from the intrinsic
spectral properties of gamma-ray sources (such as active galactic nuclei, AGN) is the primary
challenge to interpret observations of these objects. We present a Bayesian and Markov Chain Monte
Carlo approach to simultaneously infer parameters characterizing the EBL and the intrinsic spectra
in a combined fit of a set of sources, which has the advantage of easily incorporating the
uncertainties of both sets of parameters into one another through marginalization of the posterior
distribution. Taking a sample of synthetic blazars observed by the ideal CTA configuration, we
study the effects on the EBL constraints of combining multiple observations and varying their
exposure. We also apply the methodology to a set of 65 gamma-ray spectra of 36 different AGNs
measured by current Imaging Atmospheric Cherenkov Telescopes, using Hamiltonian Monte Carlo as a
solution to the difficult task of sampling in spaces with a high number of parameters. We find
robust constraints in the mid-IR region while simultaneously obtaining intrinsic spectral
parameters for all of these objects. In particular, we identify Markarian 501 (Mkn 501) flare data
(HEGRA/1997) as essential for constraining the EBL above 30 μm.