A path integration procedure based on Gauss–Legendre integration scheme is developed to analyze probabilistic solution of nonlinear vibration energy harvesters (VEHs) in this paper. First, traditional energy harvesters are briefly introduced, and their nondimensional governing and moment equations are given. These moment equations can be solved through the Runge–Kutta and Gaussian closure method. Then, the path integration method is extended to three-dimensional situation, solving the probability density function (PDF) of VEH. Three illustrative examples are considered to evaluate the effectiveness of this method. The effectiveness of nonlinearity of traditional monostable VEH is studied. The bistable VEH is further studied too. At the same time, equivalent linearization method (EQL) and Monte Carlo simulation (MCS) are employed. The results indicate that three-dimensional path integration method can give satisfactory results for the global PDF, especially when solving bistable VEH problems. The results of this method have better consistency with the simulation results than those of EQL. In addition, different degrees of hardening and softening behaviors of PDFs occur when the magnitude of nonlinearity coefficient increases or the bistable VEH is considered.