Cell-free massive MIMO is a promising technology for beyond-5G networks. Through the deployment of many cooperating access points (AP), the technology can significantly enhance user coverage and spectral efficiency compared to traditional cellular systems. Since the APs are distributed over a large area, the level of favorable propagation in cell-free massive MIMO is less than the one in colocated massive MIMO. As a result, the current linear processing schemes are not close to the optimal ones when the number of AP antennas is not very large. The aim of this paper is to develop nonlinear variational Bayes (VB) methods for data detection in cell-free massive MIMO systems. Contrary to existing work in the literature, which only attained point estimates of the transmit data symbols, the proposed methods aim to obtain the posterior distribution and the Bayes estimate of the data symbols. We develop the VB methods accordingly to the levels of cooperation among the APs. Simulation results show significant performance advantages of the developed VB methods over the linear processing techniques.