Abstract: The large-scale integration of renewable energy poses challenges to the power systems, such as the shortage of flexibility resources. However, renewable energy and power electronics connected to the distribution networks can provide flexibilities. This paper proposes a novel method to determine the stochastic feasible operation region of distribution networks, considering flexibility-providing units, network operational constraints, and uncertainties in renewable energy generation and loads. Firstly, an improved method for calculating the deterministic flexibility operation region is introduced. And then, a scenario-based approach is used to model the uncertainties of loads and distributed generations. Finally, a case study demonstrates the accuracy and effectiveness of the proposed method. The case study results show that the proposed method can provide both the probability distribution of the feasible operation region and the operation regions at different confidence levels, and achieves higher accuracy with the same computational cost and avoids overly optimistic results.