We use a modified version of the Peak-Patch excursion set formalism to compute the mass and size distribution of QCD axion miniclusters from a fully non-Gaussian initial density field obtained from numerical simulations of axion string decay . We find strong agreement with N-Body simulations at significantly lower computational cost . We employ a spherical collapse model , and provide fitting functions for the modified barrier in the radiation era . The halo mass function at z = 99 has a power-law distribution M ^ { -0.68 } for small masses and M ^ { -0.35 } for large masses within the range 10 ^ { -15 } \lesssim M \lesssim 10 ^ { -9 } M _ { \odot } , with all masses scaling as ( m _ { a } / 50 \mu \mathrm { eV } ) ^ { -0.5 } . We construct merger trees to estimate the collapse redshift and concentration mass relation , C ( M ) , which is well described using analytical results from the initial power spectrum and linear growth . Using the calibrated analytic results to extrapolate to z = 0 , our method predicts a mean concentration C \sim \mathcal { O } ( \text { few } ) \times 10 ^ { 4 } . The low computational cost of our method makes future investigation of the statistics of rare , dense miniclusters easy to achieve .