In this article, a new reliability analysis algorithm is proposed to calculate the probability density function of the bearing capacity of the foundation resting on rock mass. Despite common approaches used by other investigators, four parameters with uncertainties have been adopted in this study as random variables, including GSI index, strength of intact rock (ci), intact rock constant (mi), and rock mass disturbance factor (D). In the extended Subset Simulation (SS) proposed in this study, the samples at the first stage are produced using the Monte Carlo Simulation (MCS), while at the next levels, a Markov chain based on the Metropolis-Hastings algorithm is applied to each subset. Finally, statistical parameters of the PDF of bearing capacity are discussed. The results obtained showed that (A) The SS method converges with a much smaller number of samples than those given by the MCS method; (B) Parameters UCS and GSI have the greatest effect on the bearing capacity; (C) As the coefficient of variation of the input variables increases, the value of the reliability index decreases and therefore the probability of system failure increases.; (D) When the negative coefficient of correlation is used, a decrease in the variation of bearing capacity is observed.