Adaptive Bayesian Interval Estimation for Rare Binomial Events: A Variance-Blending Calibration Framework

1. Introduction

Rare event phenomena such as emerging disease outbreaks, earthquakes, and industrial defects are typically characterised by extremely low probability of occurrence. Estimation of binomial proportions for such events poses significant challenges, particularly due to the inadequacy of conventional confidence interval methods. Classical approaches, such as the Wald interval, frequently perform poorly for these events due to their reliance on asymptotic normality, which fails as the true proportion approaches the extreme boundaries of the binomial proportions (Wilson, 1927; Agresti and Coull, 1998). The Clopper-Pearson interval also tends to be overly conservative, thus producing wide intervals that result in statistical inefficiency (Clopper and Pearson, 1934). Conversely, asymptotic approximation methods such as the Wald interval may yield substantial under coverage probabilities when the event probability is small (Wilson, 1927).

Addressing the weaknesses of the traditional confidence intervals estimations under rare-event conditions, many frequentists alternative approaches (Newcombe 1998; Agresti and Gottard, 2007; Brown et al. 2001; Krishnamoorthy et al. 2017) and Bayesian alternatives (Lee 2006; Liu, Shao, and Yang 2025) have been proposed. Although these methods have demonstrated improvements in coverage probability, especially in small-samples or extreme and rare events, they are built under the assumption of fixed error margin, which are not plausible for events with proportions at the extreme ends of the probability spectrum (0 and 1). This is because such margins fail to scale with the magnitude of the proportion being estimated, which may lead to under- or overestimated uncertainties, depending on the location of the proportion within the parameter space (Blaker 2000). This results in inaccurate estimates because the intervals are either too wide to be practically useful or too narrow to be valid, depending on the sample size and the rarity of the event (Owen and Burke, 2024).

The incorporation of adaptive margins that scale with the magnitude of the proportion into existing classical methods has greatly improved the accuracy of both coverage probability and interval width of rare events (Owen and Burke 2024). Although alternative adaptive Bayesian methods such as the adaptive Huber’s ε-contamination model that assumes that the error margin (ε) is unknown and then adaptively estimate it has been shown to perform well in estimating probabilities on a continuous scale, but it tends to underperform in small-sample or rare-events (Luo and Gao 2024). Bayesian adaptive priors’ methods have also been proposed and shown to perform well in rare events, but require the correct scale parameter tuning, which may increase the credible width when incorrectly tuned or specified (Lyles et al. 2019). As an alternative, we introduce an adaptive Bayesian variance-blending calibration framework (hereafter known as the blended approach or method).

The proposed method offers several key advantages. It significantly improves coverage accuracy by integrating information from the beta distribution via Jeffrey’s prior, the Wilson score, and a credible level tuning parameter as well as a gamma regularization parameter which enables it to maintain the targeted nominal coverage more efficiently across a wide range of sample sizes and true proportions. The use of the beta distribution and logit transformations stabilizes estimates near the extremes; thus, our method is robust for small sample sizes and extreme proportions near 0 or 1 which provides more accurate and stable intervals. Additionally, by combining multiple sources of variability, the approach avoids overly conservative or overly optimistic interval widths, resulting in intervals that are reasonably narrow without sacrificing coverage. The method therefore strikes an optimal balance between precision and coverage. The incorporation of the gamma and the credible level parameters allows the intervals to adapt to the sample size and true proportions, which accounts for sampling uncertainty and potential overdispersion. Computationally, the proposed method is fully vectorizable and memory-efficient, thus enabling seamless integration into automated workflows and statistical software. Furthermore, its modular structure allows for rapid diagnostics and tuning, facilitating iterative refinement without computational bottlenecks.

The paper makes contributions to statistical theory of estimation by introducing a blended confidence interval method that inherits asymptotically desirable properties (consistency, efficiency, and asymptotic normality) of frequentist and model-based inferences. Specifically, as n→∞, the proposed interval exhibits consistency, asymptotic normality, and nominal coverage probability. The blending mechanism is theoretically grounded in logit-scale transformations and adaptive variance blending, ensuring that the method remains robust for small samples and across extreme events. Furthermore, the paper makes practical contributions in that it is robust and adaptive to extreme and rare events, as well as a small sample size to yield high coverage with competitively narrower width, thus making it more adaptable to problems in various fields of applications.

The rest of the paper is outlined below. The empirical literature is reviewed in section Two, while section three presents the theoretical framework. In section four, the estimation/simulation procedure and practical applications of the proposed method using COVID-19 data are presented, while the results and discussions follow in section five. The paper is then summarised in section six.

2. Review of the Literature

Interval estimation techniques for proportions have evolved significantly from simple frequentist methods to more robust model-based approaches. Historically, classical methods relied on point estimates using sample proportions and large-sample confidence intervals based on the normal approximation. However, these early methods often performed poorly, especially in small sample sizes or extreme proportions near 0 or 1. Improved frequentist methods such as the Wilson score, Agresti-Coull, and exact or Clopper-Pearson intervals were developed to offer better coverage properties. Bayesian methods have also gained prominence, allowing for the incorporation of prior knowledge which thus yielding credible intervals with intuitive probabilistic interpretations. These techniques have facilitated more flexible and accurate estimation of proportions across various fields such as epidemiology, polling, and machine learning. In this section, we review some of these methodologies.

Wilson (1927) introduced the score-based confidence interval for binomial proportions by inverting the score test which offers a more accurate alternative to the Wald interval. Designed to account for binomial asymmetry, the Wilson’s interval adjusts both the centre and width of the interval, yielding a closed-form expression that improves coverage and avoids zero-width bounds when extreme proportions are observed. It performs well across moderate sample sizes and proportions but may yield under coverage in small samples or rare event, particularly when no successes are observed. The absence of a continuity correction in its original form can further affect accuracy in discrete settings. These limitations have led to the adoption of more conservative methods, such as the exact and Bayesian intervals, for improved efficiency in rare events (Brown et.al 2001).

Clopper and Pearson (1934) introduced the exact binomial confidence interval which is grounded in the inversion of cumulative binomial tests using the equal-tail rule. The approach was designed to ensure that the true binomial proportion is contained within the interval with at least the nominal confidence level for all possible values of proportions and sample sizes. This method leverages the relationship between the binomial and beta distributions to derive interval bounds that are mathematically rigorous and distributionally valid, even in small samples. However, subsequent evaluations have shown that while the Clopper-Pearson interval maintains its nominal coverage, it is often overly conservative, particularly when the observed proportion is at the extreme ends of the probability scale or when the sample size is small. This conservatism results in wider intervals than necessary, which can reduce statistical power and limit practical applications (Brown et al. 2001).

Agresti and Coull (1998) introduced a practical and robust alternative to the traditional Wald confidence interval for estimating binomial proportions, known as the Agresti-Coull interval. Their method builds upon the Wilson score interval by adding two successes and two failures to the observed data, effectively increasing the sample size by four. This adjustment stabilizes the standard error, particularly in small samples or when the true proportion is near 0 or 1. Through extensive simulation studies, the results indicated that the proposed method achieves coverage probabilities much closer to the nominal level than both the Wald and Clopper-Pearson intervals, which tend to be either too narrow or overly conservative, respectively.

The foundational works of Clopper and Pearson (1934) & Agresti and Coull (1998) highlight the trade-off between interval accuracy and conservativeness, especially in rare event scenarios, which emphasizes the need for methods that are sensitive to rare events. In this pursuit, Blaker (2000) proposed confidence curves for discrete distributions, which offer a flexible and visually intuitive way to view coverage properties across different confidence levels. Unlike traditional fixed-level intervals that provide a single confidence range, confidence curves display the plausible parameter values across a continuum of confidence levels. One of the key strengths of the method is its ability to produce shorter and less conservative intervals than the classical Clopper-Pearson approach, while maintaining nominal coverage. This is particularly advantageous in small-sample sizes, where traditional methods often yield overly wide intervals that reduce precision. However, the algorithm is computationally intensive and can produce non-contiguous intervals in very small samples (it may exclude a narrow band of values in the middle), and the intervals may still be wide when data is extremely sparse or inherently uncertain.

Brown et al. (2001) conducted a comprehensive theoretical and empirical comparison of several binomial interval estimators. They argued that the Wald interval has an erratic coverage behaviour in small and even in large samples and then propose the use alternative methods such as the Wilson score interval, Jeffreys Bayesian interval, and Agresti-Coull interval. Using Edgeworth expansions and extensive simulations, the Wald interval was shown to significantly underperform, especially when the proportion is near 0 or 1. However, the Wilson interval was shown to perform well for small samples, while the Agresti-Coull method performed better for larger samples. Therefore, the authors recommended using the Wilson interval for small samples and the Agresti-Coull method for larger ones. Although the study provides a rigorous theoretical analysis for practical guidance, it did not address interval performance under rare-event conditions.

Seung-Chun (2006) revisited the problem of binomial proportion estimation by proposing a Bayesian interval estimator based on the weighted Polya posterior. This method blends Bayesian reasoning with frequentist properties, offering a practical alternative to traditional intervals like Wald and Agresti-Coull. The weighted Polya posterior approach yields intervals that are essentially equivalent to Agresti-Coull intervals for large samples but with improved behaviour in small-sample and skewed-distribution. Seung-Chun showed that this method maintains good coverage properties even when the sample size is small or the proportion is near 0 or 1; a condition under which many standard intervals break down. The study also discussed the link between the Polya posterior and noninformative Bayesian priors, emphasising its flexibility and ease of implementation. The findings suggest that the weighted Polya posterior is a reliable and computationally efficient method for binomial interval estimation, especially in practical applications where both robustness and simplicity are valued. However, the computational complexity and sensitivity of the method to prior specification may limit its accessibility for routine use, especially in applied settings where objective priors are preferred. The method also was applied to fixed margins, thus the need to assess how it performs under adaptive margins, which may provide opportunity for improving the accuracy of the estimates.

Ogura and Yanagimoto (2018) refined the Bayesian credible interval estimation for binomial proportions, targeting scenarios where the true proportion is small. The method’s innovation lies in the application of a logit transformation to the binomial proportion before constructing the highest posterior density (HPD) intervals. Using the transformed logit-scale, the method stabilises posterior densities and reduces the interval asymmetry which is commonly observed near the boundaries of the binomial domain. Simulation results demonstrate that the proposed transformation-based intervals yielded improved coverage and narrower width compared to traditional beta-based Bayesian intervals, especially when estimating proportions below 0.05 with modest sample sizes. The authors recommend this approach for rare-event analysis, since the logit transformation improves interpretability and numerical stability. However, the method uses a fixed-margin framework, with coverage and width derived from a static posterior formulation. It also assumes a continuous approximation and may complicate direct interpretation.

Castro et al. (2019) conducted a comparative study of interval estimation methods for binomial proportions using both Bayesian and frequentist frameworks. They evaluated the performance of equal-tailed and highest posterior density (HPD), Bayesian intervals constructed using Beta priors, alongside frequentist intervals such as Wald, Wilson, and Clopper–Pearson. Using extensive simulations across varying sample sizes and true proportions, the authors evaluated each method based on empirical coverage and expected interval width. Their results showed that Bayesian intervals-particularly HPD intervals, generally maintain better coverage stability and avoid boundary problems. In contrast, the Wald interval consistently underperforms, whereas Clopper–Pearson is too conservative. The study recommended Bayesian HPD intervals for small to moderate sample sizes and scenarios with skewed proportions. A notable strength of the method is its balanced comparison across events and samples; however, the study was limited by its reliance on noninformative priors and does not explore adaptive prior tuning or empirical Bayes approaches. This restricts the exploration of more data-responsive Bayesian intervals.

Liu et al. (2025) investigated interval estimation methods for extreme event proportions, where the binomial success probability is very low. The authors evaluate exact, approximate, and Bayesian intervals using both simulated and real-world data, including applications in environmental and reliability studies. Bayesian intervals, particularly those using noninformative or weakly informative priors, consistently outperform others in balancing coverage accuracy and interval width. The authors recommend Bayesian approaches for rare-event inference because of their flexibility and better performance. However, the study noted that theoretical justification for some Bayesian intervals is limited and that the results may be sensitive to prior choice, which could affect robustness in practice.

In contemporary times, efforts to improve estimates for coverage probabilities and interval widths have led to the shift from fixed error margins as used in traditional Bayesian and classical methods to adaptive error margins. For example, Owen and Burke (2024) focused on interval estimation for rare-event binomial proportions, where traditional methods often fail to provide more accurate estimates because of their overly wide intervals or failure to meet nominal coverages. To address these shortcomings, the authors introduced the concept of relative margin of error, which scales the interval width with the magnitude of the rare event proportion. This adjustment ensures that the interval remains accurate even when the proportion is of the order of 0.000006 to 0.1. Through simulations and analytical derivations, they showed that relative error-based criteria offer more interpretable and consistent intervals for rare events. The Wilson and Agresti-Coull intervals with adaptive margins performed better than Wald and Clopper-Pearson intervals, especially when the true proportion is close to zero. The authors recommended adopting relative margin of errors in rare event to avoid misleading conclusions. However, the study is limited to only four methods and does not explore Bayesian alternatives, nor does it address performance in moderate or large proportions, which may restrict its generalizability.

In the same year, Luo and Gao (2024) also theoretically investigated the construction of confidence intervals that are simultaneously robust and adaptive under Huber’s ε-contamination model, which allows for a small proportion of observations to deviate arbitrarily from an assumed parametric distribution. The authors established nontrivial lower bounds indicating that it is statistically infeasible to construct confidence intervals that adapt to unknown contamination levels without incurring significant inflation in interval width. Specifically, they showed that any attempt to achieve adaptive robustness-where inference remains valid across a continuum of possible model deviations-must pay an exponential price in the form of reduced efficiency. To address these limitations, the authors proposed a novel approach that relies on uniform control of uncertainty across all quantiles of the empirical distribution. Their method yielded confidence intervals with minimax-optimal coverage and near-optimal width, even when the contamination magnitude is unspecified. The construction is generalizable beyond the Gaussian location model and extends to other parametric location families and robust hypothesis testing contexts. Although the method is theoretically optimal under contamination, its dependence on distributional shape limits its application in discrete settings like binomial proportions, where contamination models are less standard. Thus, although powerful in continuous models, the approach may be less practical for binomial inference, especially in small-sample or rare-event contexts.

In the spirit of adaptive corrections of margins, Lyles et al. (2019) also explored calibrated Bayesian credible intervals methods for binomial proportions using Beta prior distributions that are adaptively tuned to control interval performance across the probability scale. The authors developed and evaluated a class of Beta priors (κ, 1−κ), where κ is a prior weight that varies depending on the target parameter, sample size and desired calibration result. Through simulation, they compare these calibrated intervals to standard Bayesian intervals and frequentist counterparts, focussing primarily on coverage probability and credible interval width. Their findings show that calibrated intervals maintain near-nominal coverage across a broad range of binomial parameters, outperforming conventional equal-tailed intervals in small sample or near-boundary probabilities. The strength of their study lies within its emphasis on empirical calibration and formal assessment of prior influence on interval behaviour. However, the method requires optimal selection of the tuning parameter, which may increase the credible width and reduce the coverage probability near the boundaries when incorrectly tunned or specified.

3. Theoretical Framework

In this section, we provide details of the traditional Bayesian interval estimation approach and then extend it to an adaptive version by incorporating our proposed margin correction and variance blending method. We also provide the methodology for the Wilson’s approach as well as methods for assessing the comparative performance of the three methods using coverage probability and confidence or credible widths. Accurate interval estimation of binomial proportions is fundamental in statistical inference. Bayesian methods offer probabilistically coherent intervals derived from the posterior distribution. When the data is binomially distributed, Bayesian updating with conjugate Beta priors produces Beta-distributed posteriors, facilitating straightforward credible interval computation. However, traditional credible intervals may be too conservative or too narrow, especially for small samples with extreme proportions. To address this, we introduce a blended variance calibration framework with tuned credible level within a Bayesian framework to enable adaptive uncertainty quantification.

3.1. Standard Bayesian Credible Interval Using the Beta Distribution and Jeffrey’s Prior

Consider n sample size (number of trials), k: number of successes, where$k\sim Binomial\left(n,p\right)$ and p, the true underlying success probability. The goal is to estimate $\left(\widehat{p}\right)$using our proposed method and the benchmark methods, the Standard Bayesian method with Jeffrey’s prior and the Wilson’s score interval.

3.1.1. Formulation

Assume that p follows the prior distribution of Jeffrey (Jeffreys 1961) such that

$$\left(p|k\right)~Beta\left({\displaystyle \frac{1}{2}}+k,{\displaystyle \frac{1}{2}}+n-k\right)$$

(1)

By denoting $a={\displaystyle \frac{1}{2}}+k$ and $b={\displaystyle \frac{1}{2}}+n-k$, and a credible level $\left(1-\alpha \right)$, then the $100\left(1-\alpha \right)\%$credible interval is defined by equation (2) as

$$\left[Bet{a}^{-1}\left({\displaystyle \frac{\alpha }{2}},a,b\right),Bet{a}^{-1}\left(1-{\displaystyle \frac{\alpha }{2}},a,b\right)\right]$$

(2)

where $Bet{a}^{-1}\left(.,a,b\right)$is the inverse CDF (percent point function) of the Beta distribution. Given the posterior mean and variance of the $Beta\left(a,b\right)$ as,

$$E\left(p\right)={\displaystyle \frac{a}{\left(a+b\right)}}\mathrm{and}Var\left(p\right)={\displaystyle \frac{ab}{{\left(a+b\right)}^2\left(a+b+1\right)}},$$

(3)

yield a credible width given by

$$Bet{a}^{-1}\left({\displaystyle \frac{\alpha }{2}},a,b\right)-Bet{a}^{-1}\left(1-{\displaystyle \frac{\alpha }{2}},a,b\right)$$

(4)

This is a quantile-based width and, although it narrows as n→∞, it does not scale directly with posterior variance. Quantiles reflect cumulative distribution, not local spread. The standard Bayesian interval shrinks with sample size but does not scale proportionally to posterior variance but reflects fixed quantile spreads instead. The sensitivity analysis in the below explains why the standard Bayesian credible interval width does not scale proportionally with the variance.

3.1.2. Non-Proportionality Scaling of the Bayesian Credible Width with Variance

Let the quantile function$Q_q\left(a,b\right)$ be implicitly defined by:

$$Q_q\left(a,b\right)={\displaystyle {\int }_0^{Q_q\left(a,b\right)}f\left(t\right)}dt,\mathrm{where}f\left(t\right)={\displaystyle \frac{t^{a-1}{\left(1-t\right)}^{b-1}}{B\left(a,b\right)}}$$

(5)

The definition in equation (5) suggests that $Q_q\left(a,b\right)$ is the value such that the cumulative distribution function (CDF) equals q i.e. $F\left(Q_q\right)=q$. Using the chain rule, we can write.

$${\displaystyle \frac{d}{dq}}F\left(Q_q\right)=\left({\displaystyle \frac{dF}{dt}}\right)\left({\displaystyle \frac{d{Q}_q}{dq}}\right)=f\left(Q_q\right){\displaystyle \frac{d{Q}_q}{dq}}$$

(6)

However, since$F\left(Q_q\right)=q$, substitution leads to ${\displaystyle \frac{d}{dq}}q=1\Rightarrow f\left(Q_q\right){\displaystyle \frac{d{Q}_q}{dq}}=1$. Therefore

$${\displaystyle \frac{d{Q}_q}{dq}}={\displaystyle \frac{1}{f\left(Q_q\right)}}$$

(7)

Substitution of equation (7) into equation (5) yields

$${\displaystyle \frac{d{Q}_q}{dq}}={\displaystyle \frac{1}{f\left(Q_q\right)}}={\displaystyle \frac{B\left(a,b\right)}{Q_q^{a-1}{\left(1-Q_q\right)}^{b-1}}}$$

(8)

Now, let $q_0$be the central quantile and ${\Delta }_q={\displaystyle \frac{\alpha }{2}}$ be an $100\left(1-\alpha \right)\%$interval, then the lower and upper quantiles are define by $Q_{\mathcal{l}}=Q_{q_0-{\Delta }_q},Q_u=Q_{q_0+{\Delta }_q}$and let the interval width be $w=Q_u-Q_{\mathcal{l}}$. Using a first-order Taylor expansion it can be shown that

$$Q_{q_0\pm {\Delta }_q}\approx Q_{q_0}\pm {\left.{\displaystyle \frac{d{Q}_q}{dq}}\right|}_{q_0}.{\Delta }_q\mathrm{and}W\approx 2{\left.{\displaystyle \frac{d{Q}_q}{dq}}\right|}_{q_0}.{\Delta }_q$$

(9)

Therefore, by subsequently substituting the quantile derivative, the interval width becomes

$$w\approx 2{\Delta }_q{\displaystyle \frac{B\left(a,b\right)}{Q_{q_0}^{a-1}{\left(1-Q_{q_0}\right)}^{b-1}}}$$

(10)

From the quantile derivative in (8)

$${\left.{\displaystyle \frac{d{Q}_q}{dq}}\right|}_{q_0}\to \infty \mathrm{if}a<0\sin \mathrm{ce}{Q}_{q_0}^{a-1}\to 0\mathrm{as}{Q}_{q_0}\to 0$$

Similarly,

$${\left.{\displaystyle \frac{d{Q}_q}{dq}}\right|}_{q_0}\to \infty \mathrm{if}b<0\sin \mathrm{ce}{\left(1-Q_{q_0}\right)}^{b-1}\to 0\mathrm{as}{Q}_{q_0}\to 1$$

The implications of these are that near the boundaries (0 and 1), the quantile function becomes extremely sensitive to small changes in q, especially when a < 1 or b < 1 leads to distorted interval widths in the tails. Therefore, the width behaves nonlinearly; more especially at the boundaries, given the quantile value $Q_{q_0}$.This explains why the intervals stretch in the tails and the width does not scale proportionally with variance.

3.2. Wilson Score Interval

Given $z={\Phi }^{-1}\left(1-{\displaystyle \frac{\alpha }{2}}\right)$, the Z-score for the confidence level $\left(1-\alpha \right)$ and $\widehat{p}={\displaystyle \frac{k}{n}}$, where k and n as defined previously, the Wilson Score interval (Wilson 1927) is formulated as

$$\left[{\displaystyle \frac{\left(\widehat{p}+\frac{z^2}{2n}\right)-z\sqrt{\frac{\widehat{p}(1-\widehat{p})}{n}+\frac{z^2}{4n^2}}}{1+\frac{z^2}{n}}},{\displaystyle \frac{\left(\widehat{p}+\frac{z^2}{2n}\right)+z\sqrt{\frac{\widehat{p}(1-\widehat{p})}{n}+\frac{z^2}{4n^2}}}{1+\frac{z^2}{n}}}\right]$$

(11)

3.3. Adaptive Bayesian-Logit-Scaled Variance-Blending Calibration Framework

Here we present a detailed formulation of our proposed Adaptive Bayesian-Logit-Scaled variance-blending calibration framework. We then prove its width’s proportionality to the margin of error and how the margin of error also adapts to the sample size by assessing its asymptotic properties.

3.3.1. Formulation

Using the same Jeffrey’s prior Beta posterior, we transform the lower and the upper limits of the confidence by first scaling it such that

$$\mathcal{l}=Bet{a}^{-1}\left({\displaystyle \frac{1-\varphi }{2}},a,b\right)\mathrm{and}u=Bet{a}^{-1}\left(1-{\displaystyle \frac{\left(1-\varphi \right)}{2}},a,b\right)$$

where $\varphi$ is the credibility multiplier, representing the coverage level and it determines how wide the interval is based on the desired confidence or credibility level. It therefore serves as a regularization parameter, and it allows calibration (via optimisation) to shrink extreme intervals to stabilize the estimates. Now let transform the interval on a logit scale as follows,

$$Logi{t}_L=\log \left({\displaystyle \frac{\mathcal{l}}{1-\mathcal{l}}}\right)\mathrm{and}Logi{t}_U=\log \left({\displaystyle \frac{u}{1-u}}\right)$$

We now compute the centre and the Half-Width as

$$c={\displaystyle \frac{Logi{t}_L+Logi{t}_U}{2}}\mathrm{and}h={\displaystyle \frac{Logi{t}_U-Logi{t}_L}{2}}$$

(12)

Definition 1. Given the sample variance$\left({\mathrm{var}}_{sample}\right)$, model-based variance$\left({\mathrm{var}}_{\gamma }\right)$ and the Wilson variance approximation$\left(MO{E}_{Wilson}\right)$ defined by

$${\mathrm{var}}_{sample}={\displaystyle \frac{1}{n\widehat{p}\left(1-\widehat{p}\right)+ϵ}},{\mathrm{var}}_{\gamma }={\left({\displaystyle \frac{\theta }{\sqrt{n}}}\right)}^2\mathrm{and}MO{E}_{Wilson}={\left(z\sqrt{{\displaystyle \frac{\widehat{p}(1-\widehat{p})}{n}}+{\displaystyle \frac{z^2}{4n^2}}}\right)}^2$$

The blended variance and the corresponding standard error are given by,

$${\mathrm{var}}_{Blend}={\mathrm{var}}_{\gamma }+{\mathrm{var}}_{sample}+MO{E}_{Wilson}\mathrm{and}S{E}_{Blend}=\sqrt{{\mathrm{var}}_{Blend}}$$

Where $ϵ=1\times {10}^{-8}$is a correction factor that prevents division by zero.

The variance blending technique is informed by frequentist variance decomposition, Bayesian regularization via prior-informed penalty terms, and decision-theoretic risk aggregation (See Blei et al. 2017; Young et al. 2005; Bayarri et al. 2004 and Wilson 1927). The sampling variance captures the uncertainty arising purely from sampling, which is derived from the Fisher Information on logit scale. It becomes large when the estimate of the proportion is near 0 or 1 or when the sample size is small. The model-based variance, also known as the adaptive Gamma regularization, acts as a stabiliser that damps erratic behaviour in small samples or when the proportion is extreme. Regularization is further enhanced by tuning the credible level parameter to balance coverage probability and interval width. The Wilson margin penalty helps stabilize intervals near the boundaries and in low-proportion regions with more conservative width.

Now, substituting the blended standard error in definition 1 into the half-width in equation (12) yields the adaptive margin of error in equation (13) below.

$$\delta =h\left({\displaystyle \frac{S{E}_{Blend}}{\max \left(S{E}_{Blend},ϵ\right)}}\right)$$

(13)

We now adjust the logit interval on the probability scale to obtain the credible interval as

$$\left[\max \left(0,{\displaystyle \frac{1}{1+e^{-\left(c-\delta \right)}}}\right),\min \left(1,{\displaystyle \frac{1}{1+e^{-\left(c+\delta \right)}}}\right)\right]$$

(14)

3.3.2. Proportionality Scaling of the Credible Width with the Variance

Let ${\Delta }_{q^{*}}={\displaystyle \frac{1-\varphi }{2}}$ , the lower bound: ${\mathcal{l}}^{*}=Q_{{\displaystyle \frac{1-\varphi }{2}}}\left(a,b\right)$ and the upper bound:$u^{*}=Q_{1-{\displaystyle \frac{\left(1-\varphi \right)}{2}}}\left(a,b\right)$, then ${\mathcal{l}}^{*}=Q_{q_0-{\Delta }_{q^{*}}},u^{*}=Q_{q_0+{\Delta }_{q^{*}}}$. The Beta quantile function $Q_{q^{*}}(a,b)$ is implicitly defined as

$$Q_{q^{*}}\text{}(a,b)={\displaystyle {\int }_0^{Q_{q^{*}}(a,b)}{f}^{*}\left(t\right)}dt,\mathrm{where}{f}^{*}\left(t\right)={\displaystyle \frac{t^{a-1}{\left(1-t\right)}^{b-1}}{B\left(a,b\right)}}$$

(15)

Applying first-order Taylor’s expansion to equation (15) yields

$$Q_{q_0\pm {\Delta }_{q*}}\approx Q_{q_0}\pm {\left.{\displaystyle \frac{d{Q}_{q*}}{d_{q*}}}\right|}_{q_0}.{\Delta }_{q*}\mathrm{and}W*\approx 2{\left.{\displaystyle \frac{d{Q}_{q*}}{d_{q*}}}\right|}_{q_0}.{\Delta }_{q*}$$

(16)

The quantile derivative function is obtained as

$${\displaystyle \frac{d{Q}_{q*}}{dq*}}={\displaystyle \frac{1}{f\left(Q_{q*}\right)}}={\displaystyle \frac{B\left(a,b\right)}{Q_{q*}^{a-1}{\left(1-Q_{q*}\right)}^{b-1}}}$$

(17)

Therefore, after the substitution of ${\Delta }_{q^{*}}={\displaystyle \frac{1-\varphi }{2}}$ , the width in equation (12) becomes,

$$W*\approx {\left.\left(1-\varphi \right){\displaystyle \frac{d{Q}_{q*}}{d_{q*}}}\right|}_{q_0}$$

(18)

Equation (18) implies that Width is proportional to $\left(1-\varphi \right)$ , thus increasing $\varphi$ reduces the interval width linearly. To further show that the width is proportional to posterior variance, we transform the interval to logit scale by letting

$$\mathcal{l}\left(q\right)=\log \left({\displaystyle \frac{Q_{q*}}{1-Q_{q*}}}\right)$$

The application of the first-order Taylor’s expansion to$\mathcal{l}\left(q\right)$ yield

$$\mathcal{l}\left(q_0\pm {\Delta }_{q*}\right)\approx \mathcal{l}\left(q_0\right)\pm {\left.{\displaystyle \frac{d\mathcal{l}}{d_{q*}}}\right|}_{q_0}.{\Delta }_{q*}$$

,

Where at $Q*=q_0$

$${\displaystyle \frac{d\mathcal{l}}{dq}}=\left({\displaystyle \frac{d\mathcal{l}}{dQ*}}\right)\left({\displaystyle \frac{dQ*}{dq*}}\right)=\left({\displaystyle \frac{1}{Q*\left(1-Q*\right)}}\right)\left({\displaystyle \frac{dQ*}{dq*}}\right)={\displaystyle \frac{B\left(a,b\right)}{Q^a*\left(1-Q^b*\right)}}$$

Therefore, by the substitution of ${\displaystyle \frac{d\mathcal{l}}{dq}}$ and ${\Delta }_{q^{*}}={\displaystyle \frac{1-\gamma }{2}}$, the width on the logit scale becomes

$$W_{Logit}^{*}\approx \left(1-\varphi \right){\displaystyle \frac{B\left(a,b\right)}{Q^a*\left(1-Q^b*\right)}}$$

(19)

By the definition of margin of error, our adjusted margin of error on the logit scale becomes

$$MO{E}_{\log it}\approx \left(1-\varphi \right){\displaystyle \frac{B\left(a,b\right)}{Q^a*\left(1-Q^b*\right)}} \sqrt{ {\displaystyle \frac{1}{n\widehat{p}\left(1-\widehat{p}\right)+ϵ}}+{\left({\displaystyle \frac{\theta }{\sqrt{n}}}\right)}^2+{\left(z\sqrt{{\displaystyle \frac{\widehat{p}(1-\widehat{p})}{n}}+{\displaystyle \frac{z^2}{4n^2}}}\right)}^2}$$

(20)

$$\Rightarrow W_{Logit}^{*}\propto \left(1-\varphi \right) \sqrt{ {\displaystyle \frac{1}{n\widehat{p}\left(1-\widehat{p}\right)+ϵ}}+{\left({\displaystyle \frac{\theta }{\sqrt{n}}}\right)}^2+ {\left(z\sqrt{{\displaystyle \frac{\widehat{p}(1-\widehat{p})}{n}}+{\displaystyle \frac{z^2}{4n^2}}}\right)}^2}$$

(21)

From equation (21), the logit scaled half-width is inversely proportional to $\varphi$ (tail probability) hence it scales proportionally to the asymptotic variance or the uncertainty.

3.3.3. Asymptotic Properties

The margin of error of Jeffrey's prior (through the Bernstein–von Mises theorem), the model-based and Wilson's method are asymptotically normal, efficient and consistent, and thus, by default, our blended method exhibits these properties. Here, we assess these properties to establish the underlying asymptotic theoretical foundation of the blended method.

(a)

Asymptotic efficiency

From definition 1, the blended variance is obtained as

$${\mathrm{var}}_{blend}=\underset{sampling}{\underbrace{{\displaystyle \frac{1}{n\widehat{p}\left(1-\widehat{p}\right)+ϵ}}}}+\underset{Gamma}{\underbrace{{\displaystyle \frac{{\theta }^2}{n}}+}}\underset{Wilson}{\underbrace{ {\left(z\sqrt{{\displaystyle \frac{\widehat{p}(1-\widehat{p})}{n}}+{\displaystyle \frac{z^2}{4n^2}}}\right)}^2}}$$

(22)

From this variance, as$n\to \infty$, each term behaves as follows: Sampling term∼ $O\left({\displaystyle \frac{1}{n}}\right)$, Gamma term ∼$O\left({\displaystyle \frac{1}{n}}\right)$ and Wilson term ∼$O\left({\displaystyle \frac{1}{n}}\right)$. Therefore, following from the Hajek–Le Cam theorem of local asymptotic normality and influence functions (Le Cam 1986), ${\mathrm{var}}_{blend}\to 0$as$n\to \infty$. Furthermore${\mathrm{var}}_{blend}\sim O\left({\displaystyle \frac{1}{n}}\right)\Rightarrow S{E}_{Blend}\sim O\left({\displaystyle \frac{1}{n}}\right)$,hence

$$\delta =h\left({\displaystyle \frac{S{E}_{Blend}}{\max \left(S{E}_{Blend},ϵ\right)}}\right)\sim O\left({\displaystyle \frac{1}{\sqrt{n}}}\right)$$

(23)

which implies that the subsequent interval width shrinks at rate $O\left({\displaystyle \frac{1}{\sqrt{n}}}\right)$confirming asymptotic efficiency.

(b)

Asymptotic consistency

By the Weak Law of Large Numbers,

$$\widehat{p}={\displaystyle \frac{k}{n}}\to p \text{be} \text{the} \text{sample} \text{proportion}.$$

Therefore, for any $\epsilon >0$,

$$\underset{n\to \infty }{\lim P}\left(\left|\widehat{p}-p\right|>\epsilon \right)=0$$

Furthermore, since we have shown that $\delta \to 0$ as$n\to \infty$and $\mathrm{logit}\left(p\right)=\log \left({\displaystyle \frac{1}{1+e^{-p}}}\right)$ is differentiable and Lipschitz on (0,1),

$$L_n,U_n\to \mathrm{logit}\left(p\right)\Rightarrow c_n\to \mathrm{logit}\left(p\right)$$

where $L_n=\max \left(0,{\displaystyle \frac{1}{1+e^{-\left(c-\delta \right)}}}\right)$, $U_n=\min \left(1,{\displaystyle \frac{1}{1+e^{-\left(c+\delta \right)}}}\right)$ and $c_n={\displaystyle \frac{Logit\left(L_n\right)+Logit\left(L_n\right)}{2}}$. Thus, the interval converges to the true population proportion as $n\to \infty$.

(c)

Asymptotic Normality

Let

$$Z_n={\displaystyle \frac{c_n-\mathrm{logit}\left(p\right)}{S{E}_{blend}}}$$

(24)

Since $c_n\to \mathrm{logit}\left(p\right)$in probability and $S{E}_{blend}\to 0$, the numerator $c_n\to \mathrm{logit}\left(p\right)$behaves like the mean of i.i.d. terms (via beta posterior), so by the Central Limit Theorem, $Z_n\stackrel{d}{\to }= N\left(0,1\right)$Hence,

$$c_n\sim N\left(\mathrm{logit}\left(p\right),{\mathrm{var}}_{blend}\right)$$

(25)

This confirms asymptotic normality in logit space. From the blended variance definition 1, it can be shown from leading order analysis that

$${\mathrm{var}}_{blend}~{\displaystyle \frac{A}{n}}\mathrm{where}A={\displaystyle \frac{1}{\widehat{p}\left(1-\widehat{p}\right)}} +{\theta }^2+z^2\widehat{p}(1-\widehat{p})$$

Therefore, following from the asymptotic normality of the observed $c_n$

$$c_n\sim N\left(\mathrm{logit}\left(p\right),{\displaystyle \frac{A}{n}}\right)$$

(26)

If we let $g\left(x\right)={\displaystyle \frac{1}{1+e^{-x}}}$ and $g^{\prime }\left(x\right)=\left({\displaystyle \frac{1}{1+e^{-x}}}\right)\left(1-{\displaystyle \frac{1}{1+e^{-x}}}\right)$, by the Delta method if $X_n\stackrel{d}{\to }N\left(\mu ,{\sigma }^2/n\right)$, then

$$g\left(X_n\right)\stackrel{d}{\to }N\left(g\left(\mu \right),{\displaystyle \frac{{\left[g\left(\mu \right)\right]}^2{\sigma }^2}{n}}\right)$$

Therefore,

$$L_n,U_n\stackrel{d}{\to }N\left(p,{\displaystyle \frac{{\left[g\left(\log it\left(p\right)\right)\right]}^{2}A}{n}}\right)$$

(27)

From the results in equation (27), the unscaled interval limits exhibit asymptotic normality, converging in distribution to normality, which implies increasing concentration around the true parameter p as the sample size grows. This classical form of asymptotic behaviour confirms the consistency of the interval bounds and validates their use as efficient estimators. Moreover, it highlights that confidence intervals become progressively narrower and more accurate in large samples, reinforcing the appropriateness of normal approximations for inference, particularly in transformed parameter spaces such as the logit scale.

4. Simulation Design

4.1. Parameters and Design

Let $\text{}k\in \{0,1,\dots ,n\}$ denote the number of observed successes in a binomial trial of size $n\in \mathbb{N}\text{}$, with estimated proportion $\widehat{p}={\displaystyle \frac{k}{n}}$, where $\widehat{p}\in \left(0,1\right)\text{}$ is the true binomial probability. For each combination of $\widehat{p}$and $n$, we perform M simulation replications following a burn-in of B samples. Each replication $i\in \left\{1,\dots ,M+B\right\}\text{}$ consists of a Bernoulli trial with Jeffreys prior-based posterior inference and three interval construction methods.

For each $i\in \left\{1,\dots ,M+B\right\}\text{}$and a combination of k and n

Step 1: Compute

$$k^{(i)}\sim Binomial(n,p)\Rightarrow {\widehat{p}}^{(i)}={\displaystyle \frac{k^{(i)}}{n}}$$

Step 2: Compute the posterior shape parameters for the Jeffreys prior$Beta(0.5,0.5)\text{}$:

$$a^{(i)}=0.5+k^{(i)},b^{(i)}=0.5+n-k^{(i)}:\text{Standard} \text{Bayesian}$$

Then, construct central$(1-\alpha )\times 100\%\text{}$posterior interval as

$${\tilde{L}}^{\left(i\right)}=Q_{{\displaystyle \frac{\alpha }{2}}}\left(a^{\left(i\right)},b^{\left(i\right)}\right),\text{}{\tilde{U}}^{\left(i\right)}=Q_{1-{\displaystyle \frac{\alpha }{2}}}\left(a^{\left(i\right)},b^{\left(i\right)}\right)$$

Step 3: Given $z={\Phi }^{-1}\left(1-{\displaystyle \frac{\alpha }{2}}\right)$, construct the Wilson interval score as follows.

Compute the centre

$${\tilde{c}}^{\left(i\right)}={\displaystyle \frac{{\widehat{p}}^{\left(i\right)}+\frac{z^2}{2n}}{1+\frac{z^2}{n}}}\mathrm{and}\mathrm{the}\mathrm{margin}{\tilde{m}}^{\left(i\right)}={\displaystyle \frac{z}{1+\frac{z^2}{n}}}\sqrt{{\displaystyle \frac{{\widehat{p}}^{\left(i\right)}(1-{\widehat{p}}^{\left(i\right)})}{n}}+{\displaystyle \frac{z^2}{4n^2}}}$$

.

Then, construct the interval

$$\left[\max \left(0,{\tilde{c}}^{\left(i\right)}-{\tilde{m}}^{\left(i\right)}\right),\min \left(1,{\tilde{c}}^{\left(i\right)}-{\tilde{m}}^{\left(i\right)}\right)\right]$$

Step 4: Construct the lower and upper bounds for the Logit margin of error as follows

Using the Jeffrey’s prior, compute the lower and upper bounds

$${\tilde{L}}^{\left(i\right)}=Q_{{\displaystyle \frac{1-\varphi }{2}}}\left(a^{\left(i\right)},b^{\left(i\right)}\right),\text{}{\tilde{U}}^{\left(i\right)}=Q_{1-{\displaystyle \frac{1-\varphi }{2}}}\left(a^{\left(i\right)},b^{\left(i\right)}\right)$$

Transform the bounds to logit scale

$${\tilde{\mathcal{l}}}_L^{\left(i\right)}=\log \left({\displaystyle \frac{{\tilde{L}}^{\left(i\right)}}{1-{\tilde{L}}^{\left(i\right)}}}\right), {\tilde{\mathcal{l}}}_U^{\left(i\right)}=\log \left({\displaystyle \frac{{\tilde{U}}^{\left(i\right)}}{1-{\tilde{U}}^{\left(i\right)}}}\right)$$

$${\tilde{c}}^{\left(i\right)}={\displaystyle \frac{{\tilde{\mathcal{l}}}_L^{\left(i\right)}+{\tilde{\mathcal{l}}}_U^{\left(i\right)}}{2}}, {\tilde{h}}^{\left(i\right)}={\displaystyle \frac{{\tilde{\mathcal{l}}}_U^{\left(i\right)}-{\tilde{\mathcal{l}}}_L^{\left(i\right)}}{2}}$$

$${\tilde{\delta }}^{\left(i\right)}=\left({\displaystyle \frac{Logi{t}_{{\tilde{U}}^{\left(i\right)}}-Logi{t}_{{\tilde{L}}^{\left(i\right)}}}{2}}\right)\left({\displaystyle \frac{S{E}_{Blend}}{\max \left(S{E}_{Blend},ϵ\right)}}\right),\tilde{S}{\tilde{E}}_{Blend}=\sqrt{{\mathrm{var}}_{Blend}}$$

$${\stackrel{⌢}{L^{\prime }}}^{\left(i\right)}=\max \left(0,{\displaystyle \frac{1}{1+e^{-\left({\tilde{c}}^{\left(i\right)}-{\tilde{\delta }}^{\left(i\right)}\right)}}}\right),{\mathrm{U}\prime }^{\left(i\right)}=\min \left(1,{\displaystyle \frac{1}{1+e^{-\left({\tilde{c}}^{\left(i\right)}-{\tilde{\delta }}^{\left(i\right)}\right)}}}\right)$$

$${L^{\prime }}^{\left(i\right)}=\max \left(0,{\displaystyle \frac{1}{1+e^{-\left(c^{\left(i\right)}-{\delta }^{\left(i\right)}\right)}}}\right),{\mathrm{U}\prime }^{\left(i\right)}=\min \left(1,{\displaystyle \frac{1}{1+e^{-\left(c^{\left(i\right)}-{\delta }^{\left(i\right)}\right)}}}\right)$$

Note that at each iteration, $\varphi$ is tunned as indicated in section 4.2.

Step 5: Computing the coverage probability and the credible width.

At each iteration M in Steps 2 to 3, calculate the coverage probability and

$$\mathrm{coverage}\mathrm{probability}={\displaystyle \frac{1}{M}}{\displaystyle \sum _{i=1}^{M}I\left(L_i\le p\le U_i\right)}$$

$$\mathrm{Width}={\displaystyle \frac{1}{M}}{\displaystyle \sum _{i=1}^M\left(U_i-L_i\right)}$$

where $I$is the indicator function defined by

$$I=\left\{\begin{array}{c}1ifp\in \left[U_i-L_i\right]\\ 0otherwise\end{array}\right.$$

and $L_i$, $U_i$ represent the corresponding lower and upper credible intervals for each of the interval estimation method.

4.2. Tuning Procedure

For each configuration $\left(n,p\right)$, the blend level ℓ is tuned using a grid search:

• Grid: $\mathcal{l}\in \left\{0.90,0.91,\dots ,0.999\right\}$For each ℓ, simulate $N=1000$ intervals.
• Compute empirical coverage and average width.
• Select ${\mathcal{l}}^{\ast }$ such that:

  $$Coverage({\mathcal{l}}^{\ast })\ge 0.95,Width({\mathcal{l}}^{\ast })=\underset{\mathcal{l}}{minWidth\left(\mathcal{l}\right)}$$

If quick validation (500-1000 trials) yields coverage < 0.95, override the tuned coverage level with ℓ=0.975. This ensures that coverages below 0.95 are minimised as much as possible.

5. Results and Discussion

The goal of this section is to estimate binomial proportion p using our proposed adaptive method to improve the coverage and precision of p for rare and extreme events that occur within the range (0.00001 and 0.99999). These ranges are typical for rare and extreme events in many sectors including healthcare, clinical trials & pharmacovigilance, psychology & behavioural economics and safety-critical engineering. The sample size range for consideration is 5 to 10000 which is typical for studies in different domains. For example, public health surveys often enrol more than 10,000 individuals to detect low-prevalence diseases.

5.1. Simulation Results

We simulate 5000 samples, where the first 1000 are use as burn-in samples to estimate proportions for rare and extreme events using our method and then compare its performance to Wilson and Jeffrey’s prior in terms of coverage probability and precision. The simulation results are summarised in Table 1, Table 2, Table 3, Table 4, Table 5 and Table 6 and Figure 1, Figure 2 and Figure 3.

5.1.1. Coverage Consistency

Across all regimes, from ultra-rare to extreme, the intervals of our blended method consistently maintain coverage probabilities close to or exceed the nominal 95% level. This is especially evident as sample size increases, where the blended method often matches or surpasses both Jeffrey’s and Wilson’s intervals. For instance, at p = 0.00001 and n = 10 000, the blended interval achieves 99.18% coverage, compared to the 90.36% under-coverages for Jeffrey and Wilson. Even in extreme tail scenarios (e.g., p = 0.99999), the blended method maintains coverage above 99% at moderate to large n, demonstrating robustness under boundary conditions. Jeffreys intervals tend to underperform at small n and extreme p, while Wilson’s intervals maintain coverage but often overcompensate with wider bounds. The results are visualized in Figure 1 below where the bleneded method tends to consistently stays above the 95% coverage probability.

5.1.2. Interval Width and Efficiency

From Table 1, Table 2, Table 3, Table 4, Table 5 and Table 6 and Figure 2, it is observed that the blended intervals tend to be wider at very small sample sizes (e.g., n = 5), which is expected due to the added gamma penalty and conservative tuning of the credible level. However, these widths remain practical and shrink rapidly as sample size increases. At n ≥ 100, the blended intervals are consistently narrower than both Jeffreys and Wilson’s intervals while preserving high coverage. For instance, at p = 0.001 and n = 1000, the blended interval width is 0.00378, compared to 0.00423 for Jeffrey and 0.00524 for Wilson. This improved efficiency reflects the adaptiveness of the blended method, which balances precision and coverage through dynamic tuning and variance blending.

5.1.3. Impact of True Probability (p)

At low true probabilities (e.g., p = 0.00001 to 0.001), the blended intervals outperform Jeffrey’s and Wilson’ intervals in both coverage and width. Jeffreys often fails to capture the true value at small n due to its symmetric prior, while Wilson’ interval maintains coverage but with inflated width. Our method, by contrast, adapts to the rarity of the event, ensuring that the interval expands appropriately to maintain coverage without excessive conservatism. As p increases to moderate levels (e.g., 0.05), all methods perform well, but the blended interval remains the most efficient. At high p (e.g., 0.99 and above), the blended method continues to excel, yielding narrower intervals and maintaining coverage even when Jeffrey’s interval begins to under-cover.

5.1.4. Sample Size Sensitivity

Smaller sample sizes (n = 5, 10) exhibit more variability and wider intervals across all methods as seen in Table 1 to 4 and Figure 2 and Figure 3. However, the blended approach maintains high coverage even under the small sample size, indicating robustness in small-sample data. For example, at p = 0.0001 and n = 10, the blended interval achieves 99.96% coverage with a width of 0.196, outperforming Jeffrey’s and Wilson’s intervals in both metrics. As sample size increases (n ≥ 100), all methods stabilize, but the blended interval consistently offers the narrowest bounds while preserving coverage. This scalability confirms that the method remains efficient across sample sizes and is particularly well-suited for large-sample studies or simulations. Notably, the Blend interval avoids the erratic discreteness-induced jumps seen in classical intervals when n is small, which is one of the strengths of the blended method.

5.1.5. Sensitivity and Robustness of the Tuning Parameter on Performance

The credible level determines the initial bounds of the Bayesian Beta interval within the blended method, before transformation and blending. Unlike fixed-level intervals, the blended approach dynamically tunes this level, typically within a range like 90% to 99.9% to achieve a target frequentist coverage (e.g., ≥95%) while minimizing interval width. This tuning is not cosmetic because it directly influences the shape, width, and efficiency of the final interval.

Coverage Calibration: The credible level is highly sensitive to the interplay between sample size and true proportion. For ultra-rare events (e.g., p = 0.00001), the method often selects a high credible level to compensate for sparse data and ensure that the interval expands enough to include the true value. Conversely, for moderate proportions and large sample sizes, the tuned level may drop to 90–93%, allowing the interval to shrink without sacrificing coverage. This sensitivity is a strength because it allows the method to self-correct based on empirical coverage feedback. If the initial tuning yields under-coverage, the method escalates the credible level (defaulting to 97.5%) to restore efficiency. This dynamic adjustment is what makes the blended method robust across different regimes where fixed-level intervals often fail.

Width Efficiency: Higher credible levels produce wider Bayesian bounds, which when transformed and blended lead to more conservative intervals. The tuning algorithm balances this trade-off by seeking the minimum credible level that still achieves the desired coverage. This ensures that the interval is not unnecessarily wide thus preserving statistical efficiency. This behaviour is evident in Table 1 to 4 and Figure 3. For example, at p = 0.00001 and n = 5, the tuned level is 98.86%, yielding a wide interval that still achieves 100% coverage. At p = 0.001 and n = 1000, the tuned level drops to 92.61%, resulting in a much narrower interval while still maintaining coverage above 97%.

Boundary Sensitivity: Near the boundaries (p ≈ 0 or p ≈ 1), the credible level tuning becomes more aggressive. This is because the logit transformation amplifies small differences near the extremes, and the method counteract this with wider initial bounds. The gamma penalty and Wilson variance blending help stabilize the interval, but the credible level remains the primary lever for controlling coverage.

The multi-panel plot in Figure 3 geometrically illustrates how the blended method balances robustness and sensitivity across extreme binomial regimes, with tuning curves that reveal deep structural responsiveness to both sample size and true proportion. For very rare and near-certain proportions, the credible levels begin conservatively high at small sample sizes, ensuring robustness against under-coverage. These curves exhibit steep initial slopes that gradually flatten, forming elongated arcs that transition from convex to nearly linear as sample size increases. This curvature reflects transformation dynamics, particularly on the logit scale, where differences near the extremes are amplified, resulting in inflated vertical displacement at low n. In contrast, mid-range proportions show gentle gradients and earlier flattening, indicating geometric stabilization and sensitivity to uncertainty. The smoothness of the curves and their alignment with the coverage probability confirm empirical responsiveness, while the mirrored trajectories across complementary proportions (e.g., p≈0.001 vs. p≈0.999) reinforce the symmetry of the transformation and principled adaptiveness. Collectively, the geometric features, steepness, curvature, flattening, and symmetry demonstrate how the method maintains high coverage while adaptively tuning credible levels across diverse regimes, independent of any fixed parameterization.

5.2. Practical Applications Using COVID-19 Data

To assess the performance of the blended method in comparison to the competing methods with applications to real life data, mortality and recovery rate for COVID-19 outbreak was used. Data from nine countries, Western Sahara, Ghana, South Africa, Australia, Brazil, Germany, France and the USA were used. These data were selected to mirror rare, extreme nature or small sample nature. The data was sourced from https://www.worldometers.info/coronavirus/. For each country, mortality and recovery rate were estimated using the blended method and the competing models. The results are reported in Table 5 to Table 6.

Across the countries, the mortality rate intervals for the blended method consistently achieved superior coverage, particularly in large samples and extremely low rates. For example, for Australia (true proportion = 0.00206), blended coverage reached 97.3%, exceeding the nominal 95% threshold and outperforming Jeffrey’s and Wilson. Similar patterns were observed for South Africa and Germany, where the blended intervals-maintained coverage above 97%. Although slightly wider in some mid-sized samples (e.g., Ghana), the blended method’s prioritization of coverage over smaller width appears justified in rare event where underestimation poses substantial epidemiological risk. In small samples such as Western Sahara (n = 10), the blended intervals preserved high coverage (98.7%) while achieving the narrowest bounds among the three methods. Jeffrey’s intervals remained stable, whereas Wilson’s intervals exhibited under coverage (92.9%), indicating sensitivity to small-sample size bias and boundary truncation.

Recovery rate estimation revealed good performance across all methods in large samples. However, blended intervals again demonstrated improved coverage exceeding 97% for Brazil and Germany. Interval widths converged across methods in large samples (width≤0.00006), yet the blended intervals retained a marginal advantage in coverage. For Ghana (true proportion = 0.90), Wilson coverage dropped to 92.7%, while blended and Jeffrey’s intervals maintained nominal levels, underscoring the blended method’s robustness near the boundary. The blended method’s adaptive tuning emerged as a key strength. This dynamic behaviour enables the method to maintain nominal coverage across heterogeneous epidemiological contexts without excessive conservatism.

The demonstrated efficiency of the blended interval method in both rare-event and boundary-extreme events supports its integration into public health surveillance systems, particularly for real-time estimation of mortality and recovery rates. Its capacity to preserve coverage in sparse data environments makes it suitable for early disease outbreak detection and low-incidence monitoring, while its stability in high-proportion contexts ensures accurate reporting of recovery rates of disease outbreaks. Adoption of the blended method in reporting frameworks could enhance the interpretability and credibility of epidemiological risks thus informing risk communication, resource allocation, and intervention prioritization. Moreover, its modular tuning architecture aligns with adaptive surveillance strategies which allows for calibration sensitivity in evolving epidemic outbreaks.

6. Conclusions

This study introduces an adaptive Bayesian variance-blending calibration framework method for estimating binomial proportions of rare and extreme events. The performance of the method in comparison with classical Jeffreys and Wilson Score intervals across different sample sizes (n=5 to 10,000) under different true proportions was assessed through extensive simulations and real data applications. The blended method was shown to be robust and adaptive across sample sizes and true probabilities. It consistently achieves or exceeds the nominal 95% coverage, particularly in small-sample, rare events and extreme events where traditional methods often struggle. For instance, at p=0.00001 and n=10,000, the blended interval attained 99.18% coverage, significantly outperforming Jeffrey’s and Wilson Score intervals. Even at small sample sizes (n=5), the method maintains high coverage and avoids the erratic behaviour typical of discrete classical intervals.

In terms of interval width, the blended method balances conservatism and efficiency through dynamic credible level tuning. While intervals are wider at small n due to conservative initialization, they shrink rapidly with increasing sample size. At moderate and large n, blended intervals are consistently narrower than Jeffreys and Wilson while preserving coverage. This adaptivity is driven by a tuning mechanism that adjusts the credible level based on the observed data and empirical coverage feedback. High credible levels are selected for ultra-rare events to ensure efficiency, while lower levels are used for moderate proportions to enhance precision. Geometric analysis of the tuning curves reveals steep, convex-to-linear transitions and mirrored symmetry across rare and extreme events spectrum which emphasizes the blended method’s adaptivity to sample sizes and events. These features reflect the method’s ability to stabilize estimates near boundaries through aggressive credible level tunning and variance blending.

Theoretically, the blended method bridges Bayesian and frequentist paradigms by using empirical feedback to calibrate credible levels to guarantee near or above nominal coverage with narrow or conservative width. Its adaptivity and boundary sensitivity offer a competitive alternative to fixed-level intervals, particularly in rare and extreme events. Practically, the method’s ability to adaptively shrink or expand based on sample size and true proportion makes it highly suitable for applications where sample size limitations and rare or extreme event proportions such as rare disease prevalence, safety-critical system reliability, and early-phase clinical trials renders conventional interval estimation inaccurate. Considering these findings, we recommend the use of the blended method for interval estimations of rare and extreme events, especially when sample sizes are limited or true proportions lie near the boundaries. Future work may extend this framework to multinomial or hierarchical settings, where adaptive tuning could further enhance inference under complex data structures.