Epidural spinal cord stimulation has been reported to partially restore volitional movement and autonomic functions after motor and sensory-complete spinal cord injury (SCI). Modern spinal cord stimulation platforms offer significant flexibility in spatial and temporal parameters of stimulation delivered. Heterogeneity in SCI and injury-related symptoms necessitate stimulation personalization to maximally restore functions. However, the large multi-dimensional stimulation space makes exhaustive tests impossible. In this paper, we present a Bayesian optimization strategy for identifying personalized optimal stimulation patterns based on the participant's expressed preference for stimulation settings. We present companion validation protocols for investigating the credibility of learned preference models. The results obtained for five participants in the E-STAND spinal cord stimulation clinical trial are reported. Personalized preference models produced by the proposed learning and optimization algorithm show that there is more similarity in optimal frequency than in pulse width across participants. Across five participants, the average model prediction accuracy is 71.5% in internal cross-validation and 65.6% in prospective validation. Statistical tests of both validation studies show that the ability of the preference models to correctly predict unseen preference data is significantly greater than chance. The personalized preference models are also shown to be significantly correlated with motor task performance across participants. We show that several aspects in participants' quality of life has been improved over the course of the trial. Overall, the results indicate that the Bayesian preference optimization algorithm could assist clinicians in the systematic programming of individualized therapeutic stimulation settings and improve the therapeutic outcomes.