The Sagdeev-Zaslavski (SZ) equation for wave turbulence is analytically derived, both in terms of a generating function and of a multipoint probability density function (PDF), for weakly interacting waves with initial random phases. When the initial amplitudes are also random, a one-point PDF equation is derived. Such analytical calculations remarkably agree with results obtained in totally different fashions. Numerical investigations of the two-dimensional nonlinear Schrödinger equation (NLSE) and of a vibrating plate prove the following: (i) Generic Hamiltonian four-wave systems rapidly attain a random distribution of phases independently of the slower dynamics of the amplitudes, vindicating the hypothesis of initially random phases. (ii) Relaxation of the Fourier amplitudes to the predicted stationary distribution (exponential) happens on a faster time scale than relaxation of the spectrum (Rayleigh-Jeans distribution). (iii) The PDF equation correctly describes dynamics under different forcings: The NLSE has an exponential PDF corresponding to a quasi-Gaussian solution, as the vibrating plates, that also shows some intermittency at very strong forcings.