An improved method for stochastic estimation of voltage sags is proposed for large distribution systems with polynomial load model. The uniform Newton-style iterative format for voltage sag magnitude is developed based on the nodal power equations associated with different faults. In the proposed method， the Levenberg-Marquardt（LM） method with variable step size is adopted to calculate the voltage sag domain and ESF indices of point of common coupling（PCC） for four faults， which are divergent by the Newton's method. It is proved on IEEE 30-bus system that load model is necessary to be involved for voltage sag assessment and the proposed method is suitable for stochastic estimation of voltage sag in large distribution systems.