As a practical answer to the real questions you're addressing, such high quantiles will generally be quite sensitive to issues with model choice (especially such things as whether you model the right censoring and how heavy the tails are in the components).
But in any case - especially when dealing with high quantiles where ordinary simulation becomes impractical - that has great value in its own right; it's an interesting question from both theoretical and practical standpoints.
A couple of other approaches to this problem are (1) using the Fast Fourier Transform and (2) direct numerical integration.
One useful reference on this topic is Luo and Shevchenko (2009)$^{[1]}$.
In it they develop an adaptive direct numerical integration approach that's faster than simulation and competitive with FFT.
The more traditional approach in actuarial work was been (3) Panjer recursion, which can be found in numerous texts. Embrechts and Frei (2009)$^{[2]}$ discuss and compare Panjer recursion and FFT. (Note that both of these techniques involve discretization of the continuous distribution.)
On the other hand, doing a very unsophisticated version of simulation, and with no effort whatever made to be efficient, generating from a compound gamma-negative binomial isn't particularly onerous. This is timing on my kids' little laptop:
system.time(replicate(100000,sum(rgamma(MASS:::rnegbin(1,4,2),5,.1))))
user system elapsed
2.82 0.00 2.84
I think 2.8 seconds to generate 100$^{\,}$000 simulations of a compound distribution on a slow little laptop really isn't bad. With some effort to be efficient (of which one might suggest many possibilities), I imagine that could be made a good deal faster.
Here's the ecdf for $10^6$ simulations (which took about 29 seconds):
$\hspace{1cm}$![enter image description here](https://i.stack.imgur.com/YKlkq.png)
We see the characteristic discrete jump at zero you expect to see with a compound distribution.
[While it should be easy to make simulation a lot faster, all three alternatives mentioned here - if carried out sensibly - should be a lot faster still.]
You should note that the actuar
package supports computation with compound distributions, and offers several methods for calculation with them.
See, for example, this vignette which discusses this facility.
[Of possibly some further passing interest, note that there is an R package for the Poisson-lognormal distribution -- poilog
; if you need that distribution at some point it may be useful.]
Added in edit:
A potential quick approximation where the gamma shape parameter isn't changing -
In the gamma case, because a convolution of gammas with constant shape parameter is another gamma, you could write down the distribution of $Y|N=n$, and then evaluate the cdf and the density at a large number of grid-values at each $n$, then simply accumulate the sum directly (rather as one would for a KDE). The direct calculation only yields a lower bound to the true quantile, but if the negative binomial is not heavy tailed it should be quite rapid.
References:
[1]: Luo, X. and Shevchenko, P.V. (2009),
"Computing Tails of Compound Distributions Using Direct Numerical Integration,"
Journal of Computational Finance, 13 (2), 73-111.
[arXiv preprint available here]
[2]: Embrechts, P., and Frei, M. (2009),
"Panjer recursion versus FFT for compound distributions,"
Mathematical Methods of Operations Research, 69:3 (July) pp 497-508.
[seems to be a pre-publication version here]
- The first case was answered in detail in this question.
- One example of the second case is shown here, where the authors apply a normal approximation to the binomial distribution used in calculations of the first case.
The third case is given by Hahn and Meeker in their handbook Statistical Intervals (2nd ed., Wiley 2017):
A two-sided $100(1-\alpha)\%$ confidence interval for $x_q$, the $q$ quantile of the normal distribution, is
$$
\left[\bar{x}-t_{(1-\alpha/2;\,n-1,\,\delta)}\frac{s}{\sqrt{n}},\;\bar{x}-t_{(\alpha/2;\,n-1,\,\delta)}\frac{s}{\sqrt{n}}\right]
$$
where $t_{(\gamma;\,n-1,\,\delta)}$ is the $\gamma$ quantile of a noncentral $t$-distribution with $n-1$ degrees of freedom and noncentrality parameter $\delta = -\sqrt{n}z_{(q)}=\sqrt{n}z_{(1-p)}$.
Here, $z_{(q)}$ denotes $\Phi^{-1}(q)$, the $q$ quantile of the standard normal distribution.
For example, let's assume we drew $n=20$ samples from a normal distribution with unknown mean and standard deviation. The sample mean was $\bar{x}=10.5$ and the sample standard deviation was $s=3.19$. Then, the two-sided $95\%$ confidence interval for the $q=0.25$ quantile $x_{0.25}$ would be given by $(6.42; 9.76)$.
Here is some R
code and a small simulation to check the coverage. By changing the parameters, you can run your own simulations:
normquantCI <- function(x, conf_level = 0.95, q = 0.5) {
x <- na.omit(x)
n <- length(x)
xbar <- mean(x)
s <- sd(x)
ncp <- -sqrt(n)*qnorm(q)
tval <- qt(c((1 + conf_level)/2, (1 - conf_level)/2), n - 1, ncp)
se <- s/sqrt(n)
xbar - tval*se
}
# Simulate the coverage
set.seed(142857)
q <- 0.25 # Quantile to calculate the CI for
conf_level <- 0.95 # Confidence level
true_mean <- 100 # The true mean of the normal distribution
true_sd <- 15 # True sd of the normal distribution
sampsi <- 20 # The sample size
trueq <- qnorm(q, true_mean, true_sd) # The true quantile
res <- replicate(1e5, {
citmp <- normquantCI(rnorm(sampsi, true_mean, true_sd), conf_level = conf_level, q = q)
ifelse(citmp[1] < trueq & citmp[2] > trueq, 1, 0)
})
sum(res)/length(res)
[1] 0.95043
Best Answer
It doesn't sound like you want the difference in quantiles: the parameter you describe is the probability of the interval $[1,2]$. More generally, you can specify any interval endpoints $[x_-, x_+]$ in advance. Given parameters $(\mu, \sigma)$ of the distribution with CDF $F_{(\mu, \sigma)}$, this probability would be written
$$\theta(\mu, \sigma) = F_{(\mu, \sigma)}(x_+) - F_{(\mu, \sigma)}(x_-).$$
As such it's just a (nice, differentiable) function of the parameters and can be addressed--at least conceptually--as would any other such function.
Specifically, let $\Lambda$ be the log likelihood (a function of $(\mu,\sigma)$) with maximum value $\Lambda_0$ and let $1-\alpha$ be the desired confidence (e.g., $\alpha=0.05$ for a 95% confidence interval). To find this confidence interval you would compute the upper $1-\alpha$ percentile $c$ of a $\chi^2(1)$ distribution (e.g., equal to 3.841 for $\alpha=0.05$) and explore the values attained by $\alpha$ within the locus of all $(\mu,\sigma)$ for which $\Lambda(\mu,\sigma) \ge 2 \Lambda_0 - c$. The range of these values forms a confidence region for $\theta$.
For example, I obtained 20 iid variates from a Normal$(1, 1/2)$ distribution (which emulates the problem situation, viewing the values as logarithms). The mean of this sample was $1.016$ and its standard deviation was $0.689$. The maximum value of twice the log likelihood equals $-40.8172$. I chose $x_- = 0$ and $x_+ = \log(2) \approx 0.693$ (which to correspond to the interval $[1,2]$ for $\exp(x)$ as in the problem statement). For this interval, the value of $\theta$ is $0.247$.
The ML estimates are the coordinates of the dot in the middle where $2 \Lambda$ is largest, because the likelihood is maximized exactly when twice its logarithm is maximized. Shaded areas show contours of $2\Lambda(\mu,\sigma)$ in the $(\mu,\sigma)$ plane in intervals of $3.841/4$, so that the region of interest (by descending values of $2\Lambda$) is comprised of the red, orange, yellow, and light yellow areas, terminated by the beginning of the blue area.
The thick lines and dashed lines, labeled 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, and 0.4, are contours of $\theta(\mu,\sigma)$. (Notice how the $0.25$ contour passes through the point of maximum likelihood. This means that $0.25$ is the ML estimate of $\theta$. Because the true value of $\theta$ is $0.247$, the ML estimate of $\theta$ is almost exactly right. This was a matter of luck, not design.)
Green dots mark where the extreme values of $\theta$ are attained: they equal $0.151$ at $\mu=1.29, \sigma=0.63$ (the right-hand dot) and they equal $0.350$ at $\mu = 0.75, \sigma=0.62$ (the left-hand dot). As the pattern of $\theta$ contours shows, all other values of $\theta$ within the region of interest lie between these two values. Therefore, a 95% CI for $\theta$ is $[0.151, 0.350]$.
This figure plots the PDFs of the distributions corresponding to the two extreme values of $\theta$. (For reference, the true PDF is shown in light dashes.) The one on the right (in blue) has the larger mean. It corresponds to the right green dot in the contour plot. Its value of $\theta$ is the blue shaded area beneath it from $0$ to $\log(2)$. The one on the left (in red) has the smaller mean. It corresponds to the left green dot in the contour plot. Its value of $\theta$ is the corresponding area beneath it. Both these PDFs are (just barely) consistent with the data: their likelihoods are within $3.841/2$ of the maximum likelihood. As you move around within the region of interest in the contour, the corresponding PDFs vary. (Indeed, some of them have more extreme means than exhibited by these two and many of them have greater standard deviations.) However, none of them has any more or any less probability between $x_-$ and $x_+$ than the two PDFs shown here.
In summary, I have described a constrained optimization problem: to find the confidence limits of $\theta$, minimize and maximize $\theta$ within the region $2 \Lambda(\mu,\sigma) \ge 2 \Lambda_0 - c$. When the log likelihood is expensive to compute, this problem can be computationally expensive to solve, but at least it is straightforward.