The size distribution of supernova remnants ( SNRs ) can help to clarify the various aspects of their evolution and interaction with the interstellar medium ( ISM ) . Since the observed samples of SNRs are a collection of objects with very different ages and origin that evolve in different conditions of the ISM , statistical Monte Carlo methods can be used to model their statistical distributions . Based on very general assumptions on the evolution , we have modeled samples of SNRs at various initial and environmental conditions , which were then compared with observed collections of SNRs . In the evolution of SNRs the pressure of the ISM is taken into account , which determines their maximum sizes and lifetimes . When comparing the modeled and observed distributions , it is very important to have homogeneous observational data free from selection effects . We found that a recently published collection of SNRs in M33 satisfies this requirement if we select the X-ray SNRs with hardness ratios in a limited range of values . An excellent agreement between distributions of this subset of SNRs and the subset of modeled SNRs was reached for a volume filling-factor of the warm phase of the ISM ( partly ionized gas with n _ { H } \sim 0.2 - 0.5 ~ { } cm ^ { -3 } ;T \sim 8000 - 10000 ~ { } K ) in M33 of \sim 90 \% . The statistical distributions constructed in this way , which reproduce practically all the statistical properties of observed SNRs , allowed us to obtain one of the important parameters of M33 : the birthrate is one SNR every { 140 } - { 150 } yr , and the total number of SNRs with a shock Mach number M _ { s } \geq 2 is larger than \sim 1000 .