The estimation of high order correlation function values is an important problem in the field of quantum computation. We show that the problem can be reduced to preparation and measurement of optical quantum states resulting after annihilation of a set number of quanta from the original beam. We apply this approach to explore various photon bunching regimes in optical states with gamma-compounded Poisson photon number statistics. We prepare and perform measurement of the thermal quantum state as well as states produced by subtracting one to ten photons from it. Maximum likelihood estimation is employed for parameter estimation. The goal of this research is the development of highly accurate procedures for generation and quality control of optical quantum states.