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 a 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^{-0.6}$ for masses within the range $10^{-15}lesssim Mlesssim 10^{-10}M_{odot}$, with all masses scaling as $(m_a/50mumathrm{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 $Csim mathcal{O}(text{few})times10^4$. The low computational cost of our method makes future investigation of the statistics of rare, dense miniclusters easy to achieve.