Most of the existing masonry bridges have been in service for a significant duration, and as a result of construction limitations, these structures often exhibit intricate geometric defects. Furthermore, under prolonged loading conditions, the rheological behavior of rock can induce deformation in masonry bridges, leading to a continuously evolving stress state. Employing an idealized model for safety assessment frequently results in an overestimation of their load-bearing capacity. To accurately evaluate the load-bearing performance and remaining service life of masonry bridges, as well as to prevent safety incidents, this study employs a parametric approach to establish a two-phase numerical model of masonry bridges. In this model, cohesive elements are introduced to simulate the bonding relationship, while the distribution pattern of geometric initial defects is determined based on the theory of conditional random fields. Additionally, the rheological behavior of rock is incorporated through a custom-written Abaqus user subroutine. Building upon this foundation, the probability distribution of the load-bearing capacity of masonry bridges is reconstructed using the maximum entropy method with fractional moment constraints. The resulting outcomes are compared and validated against those obtained using the decomposition conditional correlation matrix. Finally, the effectiveness and applicability of the proposed method are demonstrated through numerical simulations and field measurements conducted on an actual bridge. The findings reveal that the method introduced in this paper adequately accounts for the stochastic nature of geometric initial defects, objectively reflects the operational performance of masonry bridges, and effectively simulates the complete failure process of such structures. Consequently, this method provides a solid basis for the safety assessment of masonry bridges.