Probabilistic load flow (PLF) calculation, as a fundamental tool to analyze transmission system behavior, has been studied for decades. Despite a variety of available methods, existing PLF approaches rarely take system control into account. However, system control, as an automatic buffer between the fluctuations in random variables and the variations in system states, has a significant impact on the final PLF result. To consider control actions' influence, this paper proposes the first analytical PLF method for the transmission grid that takes into account primary and secondary frequency controls. This method is based on a high-precision linear power flow model, whose precision is even further improved in this paper by an original correction approach. This paper also proves that if the joint probability distribution (JPD) of random variables is expressed by a Gaussian mixture model (GMM), then the JPD of system states (e.g., nodal voltages) is an infinite GMM. By leveraging this proposition, the proposed method can generate the joint PLF of the whole system, is applicable to random variables obeying any distributions, and is capable of capturing their correlation. The high accuracy and satisfactory efficiency of this method are verified on test cases scaling from 14 to 1354 buses.