We consider the problem of peak-to-average power ratio (PAPR) reduction in orthogonal frequency-division multiplexing (OFDM) based massive multiple-input multiple-output (MIMO) downlink systems. Specifically, given a set of symbol vectors to be transmitted to K users, the problem is to find an OFDM-modulated signal that has a low PAPR and meanwhile enables multiuser interference (MUI) cancellation. Unlike previous works that tackled the problem using convex optimization, we take a Bayesian approach and develop an efficient PAPR reduction method by exploiting the redundant degrees-of-freedom of the transmit array. The sought-after signal is treated as a random vector with a hierarchical truncated Gaussian mixture prior, which has the potential to encourage a low PAPR signal with most of its samples concentrated on the boundaries. A variational expectation-maximization (EM) strategy is developed to obtain estimates of the hyperparameters associated with the prior model, along with the signal. In addition, the generalized approximate message passing (GAMP) is embedded into the variational EM framework, which results in a significant reduction in computational complexity of the proposed algorithm. Simulation results show our proposed algorithm achieves a substantial performance improvement over existing methods in terms of both the PAPR reduction and computational complexity.