Residual production cross sections in spallation reactions are key data for nuclear physics and related applications. Due to the complexity of spallation reactions, i.e., a wide range of incident energies and abundant fragments being involved, it is a challenge in theory to obtain accurately and completely energy-dependent residual cross sections. With the guidance of a simplified EPAX formulas (sEPAX), the Bayesian neural network (BNN) technique is applied to form a new machine learning model (BNN + sEPAX) for predicting fragments cross sections in proton induced nuclear spallation reactions. Three types of sample dataset for measured residual production cross sections in proton-induced nuclear spallation reactions are made, i.e., D1 consists of isotopic cross sections in reactions below 1 GeV/u, D2 consists of fragments excitation functions of reactions up to 2.6 GeV/u, and D3 the hybrid of D1 and D2. With the constructed BNN and BNN + sEPAX models, the isotopic and mass cross section distributions are compared for the 356 MeV/u $^{40}$Ca + $p$ and 1 GeV/u $^{136}$Xe + $p$ reactions, and fragments excitation functions in $^{40}$Ca + $p$, $^{56}$Fe + $p$, $^{138}$Ba + $p$ and $^{197}$Au + $p$ reactions. It is found that the BNN model needs sufficient information to achieve good extrapolations, while the BNN + sEPAX model performs better extrapolations based on less information due to the physical guidance of sEPAX formulas. It is suggested that the BNN + sEPAX model provides a new approach to predict the energy-dependent residual production cross sections produced in proton-induced nuclear spallation reactions of incident energies from tens of MeV/u up to several GeV/u.