In the course of the digitization of modern production systems, a reliable parameterization of the digital twin of machining processes is essential. For example, the digital representation of milling operations enables the process parameter selection without time-consuming and expensive test series by using stability lobe diagrams (SLD). However, the parameterization of the underlying process force model with very few cutting force experiments can prevent a reliable process design, as errors in the parameterization process are propagated to the stability analysis. Therefore, a novel two-step methodology is proposed to provide probabilistic credible intervals for conventional stability lobe diagrams: First, the unknown parameters of the process force model are estimated using a Bayesian regression method. Secondly, the estimated probability distributions of the process force parameters are used to quantify the uncertainty of the stability boundary using a mechanistic process force model. The proposed methodology is particularly characterized by its low computational cost, since time-consuming and computationally expensive Monte Carlo procedures are avoided. Instead, the methodology relies on the analytical derivation of the model parameters’ posterior probability distribution and on polynomial chaos expansion (PCE) algorithms to quantify the uncertainty in the final stability analysis.