Sub-terahertz (subTHz) antennas will play an important role in the next generations of wireless communications. Since the physical antenna size is proportional to the wavelength at its targeted central operation frequency, when comes to the subTHz frequency spectrum, the antenna fabrication tolerance needs to be accurately considered during the design stage. The classic approach to studying the average performance of the design considering fabrication tolerances is through the use of the Monte Carlo (MC) Method. In this paper, we propose an adaptive polynomial chaos expansion (PCE) method for the uncertainty quantification of subTHz horn antennas with flat-top radiation patterns. The proposed method builds a surrogate model of the antenna's response to electromagnetic excitation and estimates its statistical moments with accuracy close to the reference MC method, but with a much smaller computational complexity of roughly two orders of magnitude. This complexity gain claim is supported by the presented numerical experiments.