An approximate analytical approach using complex fractional moments (CFMs) is developed for determining the transient response probability density of nonlinear oscillators with fractional derivative elements under Gaussian white noise. The CFM is defined as the Mellin transform of a probability density function. The system response is represented in the form of a pseudo-harmonic oscillation with amplitude and phase slowly varying with time. Equivalent linearization is first implemented to obtain an equivalent natural frequency and an equivalent damping coefficient of the oscillator. In this regard, simple formulas for determining them are proposed, in which both the nonlinear elements and the elastic and viscous contributions of the fractional elements are taken into account. Next, applying stochastic averaging, the response amplitude is described by a one-dimensional stochastic differential equation, and the corresponding Fokker-Planck equation is derived. The Mellin transform then converts the Fokker-Planck equation into coupled linear ordinary differential equations governing the evolution of amplitude CFMs, which are solved with a constraint corresponding to the normalization condition for a probability density. Finally, the inverse Mellin transform of the CFMs yields the amplitude probability density. The joint probability density of displacement and velocity is also obtained from the amplitude probability density in conjunction with the Jacobian of the pseudo-harmonic variable transformation for the response. Three linear and nonlinear oscillators with fractional derivatives are considered in numerical examples. For all cases, the analytical results obtained by the proposed method are in good agreement with the results of the pertinent Monte Carlo simulation.