One of the fundamental tasks in performing robust thermoacoustic design of gas turbine combustors is calculating the modal instability risk, i.e., the probability that a thermoacoustic mode is unstable, given various sources of uncertainty (e.g., operation or boundary conditions). To alleviate the high computational cost associated with conventional Monte Carlo simulation, surrogate modeling techniques are usually employed. Unfortunately, in practice, it is not uncommon that only a small number of training samples can be afforded for surrogate model training. As a result, epistemic uncertainty may be introduced by such an “inaccurate” model, provoking a variation of modal instability risk calculation. In the current study, using Gaussian process (GP) as the surrogate model, we address the following two questions: First, how to quantify the variation of modal instability risk induced by the epistemic surrogate model uncertainty? Second, how to reduce the variation of risk calculation given a limited computational budget for the surrogate model training? For the first question, we leverage on the Bayesian characteristic of the GP model and perform correlated sampling of the GP predictions at different inputs to quantify the uncertainty of risk calculation. We show how this uncertainty shrinks when more training samples are available. For the second question, we adopt an active learning strategy to intelligently allocate training samples such that the trained GP model is highly accurate particularly in the vicinity of the zero growth rate contour. As a result, a more accurate and robust modal instability risk calculation is obtained without increasing the computational cost of surrogate model training.