The nucleation of protein aggregates and their growth are important in determining the structure of the cell's membraneless organelles as well as the pathogenesis of many diseases. The large number of molecular types of such aggregates along with the intrinsically stochastic nature of aggregation challenges our theoretical and computational abilities. Kinetic Monte Carlo simulation using the Gillespie algorithm is a powerful tool for modeling stochastic kinetics, but it is computationally demanding when a large number of diverse species is involved. To explore the mechanisms and statistics of aggregation more efficiently, we introduce a new approach to model stochastic aggregation kinetics which introduces noise into already statistically averaged equations obtained using mathematical moment closure schemes. Stochastic moment equations summarize succinctly the dynamics of the large diversity of species with different molecularity involved in aggregation but still take into account the stochastic fluctuations that accompany not only primary and secondary nucleation but also aggregate elongation, dissociation, and fragmentation. This method of "second stochasticization" works well where the fluctuations are modest in magnitude as is often encountered in vivo where the number of protein copies in some computations can be in the hundreds to thousands. Simulations using second stochasticization reveal a scaling law that correlates the size of the fluctuations in aggregate size and number with the total number of monomers. This scaling law is confirmed using experimental data. We believe second stochasticization schemes will prove valuable for bridging the gap between in vivo cell biology and detailed modeling. (The code is released on https://github.com/MYTLab/stoch-agg.).