In view of the lack of an explicit expression for the stationary response probability density of generalized nonlinear systems subjected to combined harmonic and Gaussian white noise excitations, a data-driven method is proposed in this paper. The approach involves constructing an expansion expression with undetermined coefficients and determining these coefficients through solving an optimal problem. Initially, leveraging the principle of maximum entropy and the Buckingham Pi theorem, the stationary probability density of the system energy is represented in exponential form. The power of the exponential function is then expanded into a combination of basis functions of Pi groups with undetermined coefficients, constructed from system and excitation parameters, along with the system energy. Subsequently, the coefficients are determined by solving an optimal problem aimed at minimizing the residual between the expression and histogram-based estimations of the probability density of the system energy from random state data. Additionally, a sparse optimization algorithm is employed and then the explicit expression for the probability density of the system energy can be identified including system and excitation parameters. Two typical nonlinear systems, namely the Duffing oscillator and Coulomb friction system, are given to illustrate the effectiveness and accuracy of the proposed data-driven method. The identified expressions cover both resonant and non-resonant cases, showcasing the versatility and applicability of the proposed approach. Furthermore, the extensionality of the expression is thoroughly examined and discussed.