With the increasing capacity of renewable energy sources, uncertainties regarding renewable energy and other dynamic loads in integrated energy systems (IESs) are increasing. Thus, it is necessary to study the probabilistic energy flow (PEF) of IESs. However, existing PEF calculation methods such as the point estimate method (PEM) are computationally inefficient when there are many random variables and estimated points; moreover, relatively large errors can occur when the estimated points are outside their limits. Hence, this paper presents a calculation method that addresses these problems. Because there are correlations among the variables, the Nataf transformation is employed to control the correlation quickly and effectively. A model for an IES that is interconnected with natural gas and electricity systems and accounts for the uncertainties of wind plants, photovoltaic power plants, and dynamic gas loads is presented. Correlations between wind plants and photovoltaic power plants are handled using the Nataf transformation. Finally, a modified PEM is developed to solve the PEF. For situations in which the estimated points exceed their boundaries, the power transformation and equal constraint transformation methods are used. The results of time-domain simulations demonstrate the effectiveness of the proposed approach.