The probability distribution functions (PDFs) for atomic, molecular, and total gas surface densities of M33 are determined at a resolution of about 50~pc over regions that share coherent morphological properties to unveil fingerprints of self-gravity across the star-forming disk. Most of the total gas PDFs from the central region to the edge of the star-forming disk are well-fitted by log-normal functions whose width decreases radially outwards. Because the HI velocity dispersion is approximately constant across the disk, the decrease of the PDF width is consistent with a lower Mach number for the turbulent ISM at large galactocentric radii where a higher fraction of HI is in the warm phase. The atomic gas is found mostly at face-on column densities below N$_{H}^{lim}$=2.5 10$^{21}$~cm$^{-2}$, with small radial variations of N$_{H}^{lim}$. The molecular gas PDFs do not show strong deviations from log-normal functions in the central region where molecular fractions are high. Here the high pressure and rate of star formation shapes the PDF as a log-normal function dispersing self-gravitating complexes with intense feedback at all column densities that are spatially resolved. Power law PDFs for the molecules are found near and above N$_H^{lim}$, in the well defined southern spiral arm and in a continuous dense filament extending at larger galactocentric radii; this is evident in cloud samples at different evolutionary stages along the star formation cycle. In the filament nearly half of the molecular gas departs from a log-normal PDF and power laws are also observed in pre-star forming molecular complexes. The slope of the power law is between -1 and -2. This slope, combined with maps showing where the different parts of the power law PDFs come from, suggest a power-law stratification of density within molecular cloud complexes, which is consistent with the dominance of self-gravity.