The coefficient of variation α is a critical parameter in the Brownian Passage Time (BPT) model, used to quantify the variability of large earthquake recurrence intervals. In this paper, a new estimation method is proposed for α based on paleoseismic sequences across multiple faults within a given tectonic region. By integrating Monte Carlo simulations with a Bayesian framework, the method assesses the probability distribution of α without assuming that the sample average recurrence interval equals the true mean μ, thereby avoiding epistemic bias. To validate the method, 1,000,000 simulations were conducted in two study areas of differing spatial scales. In the Western Qilian Mountains-Hexi Corridor, the posterior mean of α is 0.36 (without dating uncertainty) and 0.34 (with uncertainty). Expanding the analysis to 29 faults across western China, the estimated α increases to 0.39 (without dating uncertainty) and 0.36 (with uncertainty), with substantially reduced uncertainty bounds. The results reveal that increasing the number of paleoseismic sequences significantly reduces the uncertainty in estimating α, while considering dating uncertainty has only a minor impact. The methodology provides a robust framework for deriving region-specific recurrence variability parameters and proves particularly valuable for tectonically active regions where individual fault records are sparse but collectively form comprehensive datasets across multiple fault systems.

This study presents an improved Nishenko–Buland (NB) model to address systematic biases in estimating the coefficient of variation for earthquake recurrence intervals based on a normalizing function ${\displaystyle \frac{\mathrm{T}}{{\mathrm{T}}_{\mathrm{a}\mathrm{v}\mathrm{e}}}}$. Through Monte Carlo simulations, we demonstrate that traditional NB methods significantly underestimate the coefficient of variation when applied to limited paleoseismic datasets, with deviations reaching between 30 and 40% for small sample sizes. We developed a linear transformation and iterative optimization approach that corrects these statistical biases by standardizing recurrence interval data from different sample sizes to conform to a common standardized distribution. Application to 26 fault segments across 15 major active faults in the Hetao graben system yields a corrected coefficient of variation of α = 0.381, representing a 24% increase over the traditional method (α₀ = 0.307). This correction demonstrates that conventional approaches systematically underestimate earthquake recurrence variability, potentially compromising seismic hazard assessments. The improved model successfully eliminates sampling bias through iterative convergence, providing more reliable parameters for probability distributions in renewal-based earthquake forecasting.