Thermodynamic phase stability of three elemental boron allotropes, i.e., α-B, β-B, and γ-B, was investigated using a Bayesian interatomic potential trained via a sparse Gaussian process (SGP). SGP potentials trained with data sets from on-the-fly active learning achieve quantum mechanical level accuracy when employed in molecular dynamics (MD) simulations to predict wide-ranging thermodynamic, structural, and vibrational properties. The simulated phase diagram (500-1400 K and 0-16 GPa) agrees with experimental measurements. The SGP-based MD simulations also successfully predicted that the B13 defect is critical in stabilizing β-B below 700 K. At higher temperatures, the entropy becomes the dominant factor, making β-B the more stable phase over α-B. This letter demonstrates that SGP potentials based on a training set consisting of defect-free-only systems could make correct predictions of defect-related phenomena in solid-state crystals, paving the path to investigate crystal phase stability and transitions.