Soil properties are spatially variable due to the natural deposition process. Because of this inherent spatial variability, a slope can actually fail along any potential slip surface. A single value of Factor of safety cannot account for this variation dominated the slope stability problem. Probabilistic analysis considering the spatial variability is a reasonable method to quantify the risk of the slope stability problem. Thus, in order to better simulate this variation, the theory of random field has been widely used in the slope stability problem. However, statistical outcomes derived from the probabilistic analysis will be influenced by how the random fields are generated and how the random field values are assigned to each potential slip surface. In order to investigate the extent of this influence, this study proposed a new probabilistic slope stability analysis method and compared it with the other two methods in terms of accuracy and efficiency. In this report, three different probabilistic slope stability analysis methods are presented. These methods combined the traditional limit equilibrium method of slices with random fields, which can account for the inherent spatial variability of soil properties. An exponential decaying function is used to describe the correlation structure of this spatial variation. This correlation structure is further expressed by the form of covariance matrix. Since the covariance matrix is a symmetric positive definite matrix, Cholesky decomposition is used to decompose it into the product of two triangular matrices. Because the triangular matrix is more computationally efficient, the two-dimensional random field is generated by multiplying a normal random number vector with the lower triangular matrix derived from Cholesky decomposition.