The quantitative risk assessment of landslides requires the spatial impact probability of their run-out processes to conduct vulnerability assessments and prioritize mitigation measures. Hence, a general framework for three-dimensional (3D) probabilistic landslide run-out hazard evaluation is presented in this paper. In the proposed framework, a continuum mechanics concept-based dynamic numerical model, called Massflow, is used to conduct run-out analyses. A multi-observation Bayesian back analysis method with Markov chain Monte Carlo simulation is used to calibrate the uncertainties of the input parameters based on a past landslide event, which integrates the prior information, observation information, and model bias factor to yield posterior distributions of the input parameters. Then, the posterior distributions are used as inputs for a future event to estimate the exceedance probabilities of run-out intensities (i.e., maximum run-out height and velocity) on a 3D terrain and produce run-out hazard maps for quantitative landslide risk assessments purposes. To improve the computational efficiency of the proposed framework without losing accuracy, a multi-response Kriging-based surrogate model with a sequential adaptive sampling strategy is developed for both the back analysis and forward prediction. Finally, two sequential slides of landslide CJ#8 located in the Heifangtai terrain, Gansu, China are employed to test the performance of the proposed framework, with its advantages and potential drawbacks discussed in detail. Results show that the proposed framework can effectively calibrate the model input parameters to evaluate the run-out hazard of similar landslides, and the predicted results can be employed in landslide risk assessments and mitigations.