<p>This paper investigates the problem of reconstructing a random source from statistical phaseless data for the two-dimensional Helmholtz equation. The major challenge of this problem is non-uniqueness, which we overcome through a reference source technique. Firstly, we introduce some artificially added point sources into the inverse random source system and derive phase retrieval (PR) formulas for the expectation and variance of the radiated fields. This paper rigorously analyzes the uniqueness and stability of the recovered statistics of the radiated fields. Afterwards, since the direct problem has a unique mild solution, by examining the expectation and variance of this solution and combined with the phase retrieval formulas, we derive the Fredholm integral equations to solve the inverse random source problem (IRSP). We prove the stability of the corresponding integral equations. To quantify the uncertainty of the random source, we utilize the Bayesian method to reconstruct the random source and establish the well-posedness of the posterior distribution. Finally, numerical experiments demonstrate the effectiveness of the proposed method and validate the theoretical results.</p>