  We use MCMC to sample the posterior distribution of the mass of the graviton
-- assumed here to be manifest through a Yukawa suppression of the Newtonian
potential -- by using INPOP planetary ephemerides. The main technical
difficulty is the lack of analytical formulation for the forward problem and
the cost in term of computation time for its numerical estimation. To overcome
these problems we approximate an interpolated likelihood for the MCMC with the
Gaussian Process Regression. We also propose a possible way to assess the
uncertainty of approximation of the likelihood by mean of some realization of
the Gaussian Process. At the end of the procedure, a 99.7% confidence level
threshold value is found at $1.01 \times 10^{-24} \; eV c^{-2}$ (resp.
$\lambda_g \geq 122.48 \times 10^{13} \; km$), representing an improvement of 1
order of magnitude relative to the previous estimation of Bernus et al. 2020.
Beyond this limit, no clear information is provided by the current state of the
planetary ephemerides.
