Abstract | Seismic fragility curves provide the probability of exceedance of a given damage state, should different levels of ground motion intensity be experienced at the site where the structure, or component, is located. Such curves are often derived via multiple nonlinear response history analyses (NLRHA) using sets of “suitable” ground motions that, in line with the best practice, should be consistent with the seismic hazard at the site. Based on the selected sets of records, one can estimate fragility functions that are often assumed to follow a lognormal distribution defined by two parameters, i.e., the logarithmic mean (µ) and the logarithmic standard deviation (β). Our focus is on estimating them using a state-of-the-art approach that involves hazard-consistent record selection via Conditional Spectrum and multiple stripe analysis. However, this approach usually requires many NLRHAs, with high computational costs, especially for the complex structural models typical of the nuclear industry. This study investigates the optimal number of ground motions and intensity levels required to keep the computational burden acceptable without compromising accuracy. To do so, we adopt a Bayesian framework with Markov chain Monte Carlo simulation and Metropolis–Hasting sampling. Our findings show that this approach effectively helps analysts best allocate computational resources while ensuring acceptable accuracy in estimating the probability of reaching or exceeding the considered damage states.