Understanding and accounting for uncertainty is one aspect of ensuring next-step tokamaks such as SPARC will robustly achieve their goals. While traditional Plasma OPerating CONtour (POPCON) analyses guide design, they often overlook the significant impact of uncertainties in scaling laws, plasma profiles, and impurity concentrations on performance predictions. This work confronts these challenges by introducing statistical POPCONs, which leverage Monte Carlo analysis to quantify the sensitivity of SPARC’s operating points (Creely et al 2020 J. Plasma Phys.86 5) to these crucial variables. For profiles, a physically motivated gradient-based functional form is introduced. We further develop a multi-fidelity Bayesian optimization workflow that effectively identifies operating points maximizing the probability of meeting performance goals, which gives a significant speed-up over brute force search methods. Our findings reveal that accounting for these uncertainties leads to an optimal operating point different from deterministic predictions, which balances H-mode access, confinement, impurity dilution, and auxiliary power.