This paper proposes a solution procedure to formulate an approximate joint probability density function (PDF) of a Duffing-type energy harvester system under Gaussian white noise. The joint PDF solution of displacement, velocity, and an electrical variable is governed by the Fokker-Planck (FP) equation. First, the FP equation is reduced to a lower-dimensional FP equation only about displacement and velocity by a state-space-split (SSS) method. The stationary joint PDF of displacement and velocity can be solved exactly. Then, the joint PDF of displacement, velocity, and the electrical variable can be approximated by the product of the obtained exact PDF and the conditional Gaussian PDF of the electrical variable. A parametric study is further conducted to show the effectiveness of the proposed solution procedure. The study considers weak nonlinearity, strong nonlinearity, high excitation level, and a bistable oscillator. Comparison with the simulated results shows that the proposed solution procedure is effective in obtaining the joint PDF of the energy harvester system in the examined examples.