Shales can form a complex fracture network during hydraulic fracturing, which greatly increases the stimulated reservoir volume (SRV) and thus significantly increases oil or gas production. It is therefore important to accurately predict the probability of formation of the hydraulic fracture network for shale gas exploration and exploitation. Conventional discriminant criteria are presented as the relationship curves of stress difference vs. intersection angle. However, these methods are inadequate for application in the field. In this study, an effective and quantitative prediction method relating to the probability of complex fracture network formation is proposed. First, a discriminant criterion of fracture network was derived. Secondly, Monte Carlo simulation was applied to calculate the probability of the formation of the complex fracture network. Then, the method was validated by applying it to individual wells of two active shale gas blocks in the Sichuan Basin, China. Results show that the probabilities of fracture network are 0.98 for well JY1 and 0.26 for well W204, which is consistent with the micro-seismic hydraulic fracturing monitoring and actual gas production. Finally, the method was further extended to apply for the regional scale of the Sichuan Basin, where the general probabilities of fracture network formation are 0.32–1 and 0.74–1 for Weiyuan and Jiaoshiba blocks, respectively. The Jiaoshiba block has, therefore, an overall higher probability for formation of fracture network than the Weiyuan block. The proposed method has the potential in further application to evaluation and prediction of hydraulic fracturing operations in shale reservoirs.