In the inverse problem of photoacoustic tomography (PAT), initial pressure distribution induced by the photoacoustic effect is estimated from a set of measured ultrasound data. In the recent decade, utilization of various deep learning approaches for the inverse problem of PAT have been proposed. However, many of these approaches do not provide information of the uncertainties of the reconstructed images. In this work, we propose a deep learning based approach for the Bayesian inverse problem of PAT based on variational autoencoders. The approach is evaluated using numerical simulations and compared against posterior distribution obtained using a conventional Bayesian image reconstruction approach. The approach is shown to provide rapid and accurate reconstructions with reliability estimates.