  The data-driven Bayesian model averaging is a rigorous statistical approach
to combining multiple models for a unified prediction. Compared with the
individual model, it provides more reliable information, especially for
problems involving apparent model dependence. In this work, within both the
non-relativistic Skyrme energy density functional and the nonlinear
relativistic mean field model, the effective proton-neutron chemical potential
difference $\Delta \mu^*_{\rm{pn}}$ of neutron-rich nuclei is found to be
strongly sensitive to the symmetry energy $E_{\rm{sym}}(\rho)$ around
$2\rho_0/3$, with $\rho_0$ being the nuclear saturation density. Given
discrepancies on the $\Delta \mu^*_{\rm{pn}}$-$E_{\rm{sym}}(2\rho_0/3)$
correlations between the two models, we carry out a Bayesian model averaging
analysis based on Gaussian process emulators to extract the symmetry energy
around $2\rho_0/3$ from the measured $\Delta \mu^*_{\rm{pn}}$ of 5 doubly magic
nuclei $^{48}$Ca, $^{68}$Ni, $^{88}$Sr, $^{132}$Sn and $^{208}$Pb.
Specifically, the $E_{\mathrm{sym}}(2\rho_0/3)$ is inferred to be
$E_{\mathrm{sym}}(2\rho_0/3) = 25.6_{-1.3}^{+1.4}\,\mathrm{MeV}$ at $1\sigma$
confidence level. The obtained constraints on the $E_{\mathrm{sym}}(\rho)$
around $2\rho_0/3$ agree well with microscopic predictions and results from
other isovector indicators.
