Seismic fragility models provide a probabilistic relation between ground-motion intensity and damage, making them a crucial component of many regional risk assessments. Estimating such models from damage data gathered after past earthquakes is challenging because of uncertainty in the ground-motion intensity the structures were subjected to. Here, we develop a Bayesian estimation procedure that performs joint inference over ground-motion intensity and fragility model parameters. When applied to simulated damage data, the proposed method can recover the data-generating fragility functions, while the traditionally used method, employing fixed, best-estimate, intensity values, fails to do so. Analyses using synthetic data with known properties show that the traditional method results in flatter fragility functions that overestimate damage probabilities for low-intensity values and underestimate probabilities for large values. Similar trends are observed when comparing both methods on real damage data. The results suggest that neglecting ground-motion uncertainty manifests in apparent dispersion in the estimated fragility functions. This undesirable feature can be mitigated through the proposed Bayesian procedure.