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=629$ has a power-law distribution ${M}^{\ensuremath{-}0.6}$ for masses within the range ${10}^{\ensuremath{-}15}\ensuremath{\lesssim}M\ensuremath{\lesssim}{10}^{\ensuremath{-}10}\text{ }\text{ }{M}_{\ensuremath{\bigodot}}$, with all masses scaling as $({m}_{a}/50\text{ }\text{ }\ensuremath{\mu}\mathrm{eV}{)}^{\ensuremath{-}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\ensuremath{\sim}\mathcal{O}(\mathrm{few})\ifmmode\times\else\texttimes\fi{}{10}^{4}$. The low computational cost of our method makes future investigation of the statistics of rare, dense miniclusters easy to achieve.