  Recovering credible cosmological parameter constraints in a weak lensing
shear analysis requires an accurate model that can be used to marginalize over
nuisance parameters describing potential sources of systematic uncertainty,
such as the uncertainties on the sample redshift distribution $n(z)$. Due to
the challenge of running Markov Chain Monte-Carlo (MCMC) in the high
dimensional parameter spaces in which the $n(z)$ uncertainties may be
parameterized, it is common practice to simplify the $n(z)$ parameterization or
combine MCMC chains that each have a fixed $n(z)$ resampled from the $n(z)$
uncertainties. In this work, we propose a statistically-principled Bayesian
resampling approach for marginalizing over the $n(z)$ uncertainty using
multiple MCMC chains. We self-consistently compare the new method to existing
ones from the literature in the context of a forecasted cosmic shear analysis
for the HSC three-year shape catalog, and find that these methods recover
similar cosmological parameter constraints, implying that using the most
computationally efficient of the approaches is appropriate. However, we find
that for datasets with the constraining power of the full HSC survey dataset
(and, by implication, those upcoming surveys with even tighter constraints),
the choice of method for marginalizing over $n(z)$ uncertainty among the
several methods from the literature may significantly impact the statistical
uncertainties on cosmological parameters, and a careful model selection is
needed to ensure credible parameter intervals.
