A fragility function quantifies the probability that a structural system exposed to a given hazard exceeds an undesirable limit state event conditioned on the occurrence of a hazard level. Multiple sources of uncertainty affect this function, including record-to-record variation, geometric and material properties, aging, modeling assumptions and errors, and even the analyzed dataset. This study presents a methodology for statistical model selection and uncertainty quantification of seismic fragility functions. The statistical models are created by implementing a hierarchical Bayesian framework with a sequential Monte Carlo technique. The most probable model is selected using Bayesian model selection. This model is validated through multiple metrics using predictive intervals and the Kolmogorov-Smirnov test. Then, the epistemic uncertainty is quantified as the variance of the area under the fragility functions. The methodology is implemented on a twenty-story steel benchmark model case study, demonstrating that the log-normal distribution yields superior performance relative to other models considered. Finally, further analysis of the case study demonstrates that the epistemic uncertainty is considerably reduced when using forty observations.