This paper presents the probabilistic analysis of landslides in spatially variable soil deposits, modeled by a stochastic framework which integrates the random field theory with generalized interpolation material point method (GIMP). Random fields are simulated using Cholesky matrix decomposition (CMD) method and Latin hypercube sampling (LHS) method, which represent material properties discretized into sets of random soil shear strength variables with statistical properties. The approach is applied to landslides in clayey deposits under undrained conditions with random fields of undrained shear strength parameters, in order to quantify the uncertainties of post-failure behavior at different scales of fluctuation (SOF) and coefficients of variation (COV). Results show that the employed approach can reliably simulate the whole landslide process and assess the uncertainties of runout motions. It is demonstrated that the natural heterogeneity of shear strength in landslides notably influences their post-failure behavior. Compared with a homogeneous landslide model which yields conservative results and underestimation of the risks, consideration of heterogeneity shows larger landslide influence zones. With SOF values increasing, the variances of influence zones also increase, and with higher values of COV, the mean values of the influence zone also increase, resulting in higher uncertainties of post-failure behavior.