In this paper, the probability density function (PDF) and the mean up-crossing rate of nonlinear ship rolling in random beam seas are investigated. The excitation of stationary random sea waves is approximated as a second-order linear filtered white noise. The Fokker–Planck–Kolmogorov (FPK) equation governing the probability density function of ship rolling is a four-dimensional linear partial differential equation with varying coefficients, and obtaining its exact solution is much more sophisticated. The exponential-polynomial closure (EPC) method is applied to solve the corresponding FPK equation of the system. In numerical examples, linear-plus-cubic damping model and linear-plus-quadratic damping model with three different sea states are further examined. Comparison with the equivalent linearisation (EQL) method and Monte Carlo simulated results show that the proposed procedure is effective to obtain a satisfactory probability density function solution, especially in the tail region.