Compact and highly performing photonic devices are characterized by non-intuitive geometries, a large number of parameters, and multiple figures of merit. Optimization and machine learning techniques have been explored to handle these complex designs, but the existing approaches often overlook stochastic quantities. As an example, random fabrication uncertainties critically determines experimental device performance. Here, we present a novel approach for the stochastic multi-objective design of photonic devices combining unsupervised dimensionality reduction and Gaussian process regression. The proposed approach allows to efficiently identify promising alternative designs and model the statistic of their response. Incorporating both deterministic and stochastic quantities into the design process enables a comprehensive analysis of the device and of the possible trade-offs between different performance metrics. As a proof-of-concept, we investigate surface gratings for fiber coupling in a silicon-on-insulator platform, considering variability in structure sizes, silicon thickness, and multi-step etch alignment. We analyze 86 alternative designs presenting comparable performance when neglecting variability, discovering on the contrary marked differences in yield and worst-case figures for both fiber coupling efficiency and back-reflections. Pareto frontiers demonstrating optimized device robustness are identified as well, offering a powerful tool for the design and optimization of photonic devices with stochastic figures of merit.