Stress-strength reliability is the probability that a random strength variable exceeds a random stress variable. This study develops frequentist and Bayesian inference for \(R=P(Y<X)\) when both variables follow two-parameter negative binomial (NB2) distributions. Reliability is evaluated using an adaptively truncated infinite sum with a deterministic error bound. Model parameters are estimated by maximum likelihood on the log scale. Frequentist inference includes the delta method, an equality-constrained likelihood-ratio procedure, and a parametric percentile bootstrap. Bayesian inference uses a multivariate Lindley approximation and adaptive random-walk Metropolis sampling under weakly informative priors. A frequentist Monte Carlo study used 10, 000 replications per configuration and 1000 bootstrap samples per replication. This produced aggregate unconditional coverage ranges of \(0.9443-0.9495\) for the delta method, \(0.9457-0.9466\) for the likelihood-ratio method, and \(0.9376-0.9437\) for the bootstrap. In a focused Bayesian study with 2000 replications per configuration, Markov chain Monte Carlo (MCMC) credible-interval coverage ranged from 0.9370 to 0.9575. The Lindley interval ranged from 0.8830 to 0.9675 and was least accurate in small-sample near-Poisson settings. Interval widths and estimation errors generally decreased with sample size. An automobile insurance application illustrates the framework’s practical interpretation and shows the importance of accounting for tied counts.
