The natural gas system and electricity system are coupled tightly by gas turbines in an integrated energy system. The uncertainties of one system will not only threaten its own safe operation but also be likely to have a significant impact on the other. Therefore, it is necessary to study the variation of state variables when random fluctuations emerge in the coupled system. In this paper, a multislack-bus model is proposed to calculate the power and gas flow in the coupled system. A unified probabilistic power and gas flow calculation, in which the cumulant method and Gram–Charlier expansion are applied, is first presented to obtain the distribution of state variables after considering the effects of uncertain factors. When the variation range of random factors is too large, a new method of piecewise linearization is put forward to achieve a better fitting precision of probability distribution. Compared to the Monte Carlo method, the proposed method can reduce computation time greatly while reaching a satisfactory accuracy. The validity of the proposed methods is verified in a coupled system that consists of a 15-node natural gas system and the IEEE case24 power system.