  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 interpreting 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$\mu$m.
