  Non-Gaussian statistics of the projected weak lensing field are powerful
estimators that can outperform the constraining power of the two-point
functions in inferring cosmological parameters. This is because these
estimators extract the non-Gaussian information contained in the small scales.
However, fully leveraging the statistical precision of such estimators is
hampered by theoretical uncertainties, such as those arising from baryonic
physics. Moreover, as non-Gaussian estimators mix different scales, there
exists no natural cut-off scale below which baryonic feedback can be completely
removed. We therefore present a Bayesian solution for accounting for baryonic
feedback uncertainty in weak lensing non-Gaussianity inference. Our solution
implements Bayesian model averaging (BMA), a statistical framework that
accounts for model uncertainty and combines the strengths of different models
to produce more robust and reliable parameter inferences. We demonstrate the
effectiveness of this approach in a Stage IV convergence peak counts analysis,
including three baryonic feedback models. We find that the resulting BMA
posterior distribution safeguards parameter inference against biases due to
baryonic feedback, and therefore provides a robust framework for obtaining
accurate cosmological constraints at Stage IV precision under model uncertainty
scenarios.
