This paper presents an empirical Bayesian method to estimate regularisation parameters in imaging inverse problems. The method calibrates regularisation parameters directly from the observed data by maximum marginal likelihood estimation, and is useful for inverse problems that are convex. A main novelty is that maximum likelihood estimation is performed efficiently by using a stochastic proximal gradient algorithm that is driven by two proximal Markov chain Monte Carlo samplers, intimately combining modern optimisation and sampling techniques. The proposed methodology is illustrated with an application to total-variation image deconvolution, where it compares favourably to alternative Bayesian and non-Bayesian approaches from the state of the art.