Context. Self-calibration methods with the CLEAN algorithm have been widely employed in very long baseline interferometry (VLBI) data processing in order to correct antenna-based amplitude and phase corruptions present in the data. However, human interaction during the conventional CLEAN self-calibration process can impose a strong effective prior, which in turn may produce artifacts within the final image and hinder the reproducibility of final results. Aims. In this work, we aim to demonstrate a combined self-calibration and imaging method for VLBI data in a Bayesian inference framework. The method corrects for amplitude and phase gains for each antenna and polarization mode by inferring the temporal correlation of the gain solutions. Methods. We use Stokes I data of M87 taken with the Very Long Baseline Array (VLBA) at43 GHz, pre-calibrated using the rPICARD CASA-based pipeline. For antenna-based gain calibration and imaging, we use the Bayesian imaging software resolve. To estimate gain and image uncertainties, we use a variational inference method. Results. We obtain a high-resolution M87 Stokes I image at 43 GHz in conjunction with antenna-based gain solutions using our Bayesian self-calibration and imaging method. The core with counter-jet structure is better resolved, and extended jet emission is better described compared to the CLEAN reconstruction. Furthermore, uncertainty estimation of the image and antenna-based gains allows us to quantify the reliability of the result. Conclusions. Our Bayesian self-calibration and imaging method is able to reconstruct robust and reproducible Stokes I images and gain solutions with uncertainty estimation by taking into account the uncertainty information in the data.