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 $Δμ^*_{\rm{pn}}$ of neutron-rich nuclei is found to be strongly sensitive to the symmetry energy $E_{\rm{sym}}(ρ)$ around $2ρ_0/3$, with $ρ_0$ being the nuclear saturation density. Given discrepancies on the $Δμ^*_{\rm{pn}}$-$E_{\rm{sym}}(2ρ_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ρ_0/3$ from the measured $Δμ^*_{\rm{pn}}$ of 5 doubly magic nuclei $^{48}$Ca, $^{68}$Ni, $^{88}$Sr, $^{132}$Sn and $^{208}$Pb. Specifically, the $E_{\mathrm{sym}}(2ρ_0/3)$ is inferred to be $E_{\mathrm{sym}}(2ρ_0/3) = 25.6_{-1.3}^{+1.4}\,\mathrm{MeV}$ at $1σ$ confidence level. The obtained constraints on the $E_{\mathrm{sym}}(ρ)$ around $2ρ_0/3$ agree well with microscopic predictions and results from other isovector indicators.