In computational fluid dynamics, standard geometries and flow conditions can result in complex local phenomena and interactions, such as shocks and flow separations, considerably challenging the computed solution's quality. For example, the efficient computation of supersonic compressors and turbines is challenging in turbomachinery because of complex inter-blade shock interactions that induce significant pressure loss. These critical flow structures generally require fine discretizations to be correctly computed. Their exact location is sensitive to boundary conditions and geometries, and it is usually impossible to anticipate their positions. Since uniform fine meshes would lead to overly costly simulations, alternative meshing strategies are in order. Adaptive mesh refinement (AMR) is classically employed to adjust the local mesh resolution automatically to the flow structures and construct meshes that achieve a prescribed precision while limiting the computational cost. The present work considers the case of random flow conditions sampling strategies (Monte Carlo method) to estimate some flow statistics. Adapting the mesh for each sampled condition is usually prohibitively costly while using the mesh adapted to a particular condition for all samples can yield large errors. Therefore, we propose to build a unique adapted mesh to minimize the average error. The iterative adaptation uses an estimation of the local mean error from a reduced sample set of conditions to mark the cells for refinement. In the present implementation, the cell's refinement consists of an isotropic insertion of nodes, relying on the utility adaptCells of Cassioppée. We first demonstrate the method on a one-dimensional parametric Burgers equation with known exact solutions; this simple problem enables the complete statistical characterization of the resulting errors and their sensitivity to the numerical parameters of the adaptation. Then, a two-dimensional supersonic scramjet inlet flow is considered. The model is the Euler equations solved with the cell-centered finite-volume method implemented in the CFD solver elsA. This case is widely used in the mesh adaptation literature as it yields complex shock structures and interactions. Results show the ability of the method to refine all regions requiring adaptation and reduce the mean error efficiently over the range of inlet Mach numbers.