#> Loading required package: ssanv
#> Loading required package: testthat
Klaschka and Reiczigel (2020) gave better algorithms for calculating the Blaker and minlike (i.e., Sterne) two-sided confidence intervals. They focus on the discontinuity points. Borrowing from those some of the main ideas in that paper, but using different notation, we create an algorithm to calculate the \(100(1-\alpha)\)% two-sided Poisson confidence intervals. Checks of this algorithm show that it gets nearly identical results as Klaschka and Reiczigel (2020) in all situations that were checked.
Let \(X \sim Poisson(\theta)\). Let a general two-sided p-value for testing \(H_0: \theta=\theta_0\) versus \(H_1: \theta \neq \theta_0\) at \(X=x\) be \(p_T(x,\theta_0)\). Let the associated two-sided \(100(1-\alpha)\)% confidence interval be \[\left( L_T(x, 1-\alpha), U_T(x, 1-\alpha) \right),\] where \[L_T(x, 1-\alpha) = \sup \left\{ \theta: \mbox{ for all } \theta_0 \leq \theta \mbox{ then } p_T(x, \theta_0) \leq \alpha \right\}\] and \[U_T(x, 1-\alpha) = \inf \left\{ \theta: \mbox{ for all } \theta_0 \geq \theta, \mbox{ then } p_T(x, \theta_0) \leq \alpha \right\}\].
Let \(f(i; \theta)=Pr[ X=i | \theta]\) be the Poisson probability mass function at \(X=i\). Then the ‘minlike’ p-value for testing \(H_0: \theta=\theta_0\) is \[p_m(x, \theta_0) = \sum_{i: f(i;\theta_0) \leq f(x;\theta_0)} f(i;\theta_0).\] We can write \(p_m\) in a format that is conducive for calculating the confidence interval. First, define \[\tilde{\theta}_{a,b}\] as the geometric mean of \(a+1, a+2,\ldots, b\), where \(a<b\) and both are non-negative integers. Then we can show that \(f(a; \tilde{\theta}_{a,b}) = f(b, \tilde{\theta}_{a,b})\).
Writing out \(f(i;\tilde{\theta}_{a,b})\) for \(i=a,b\), and setting them equal gives,
\[\frac{\tilde{\theta}_{a,b}^{a} e^{-\tilde{\theta}_{a,b}} }{a!} = \frac{\tilde{\theta}_{a,b}^{b} e^{-\tilde{\theta}_{a,b}} }{b!},\] and solving for \(\tilde{\theta}_{a,b}\) gives the definition of that geometric mean, \[\tilde{\theta}_{a,b} = \exp \left( \frac{1}{b-a} \sum_{i=a+1}^{a+(b-a)} \log(i) \right).\] Thus, for \(x \geq 2\) and \(2 \leq j \leq x\), we have that for \(\tilde{\theta}_{x-j,x} \leq \theta_0 < \tilde{\theta}_{x-j+1,x}\) we have \[f(x-j;\theta_0) \leq f(x; \theta_0) < f(x-j+1;\theta_0).\] So that for \(\tilde{\theta}_{x-j,x} \leq \theta_0 < \tilde{\theta}_{x-j+1,x}\), we have \[p_m(x;\theta_0) = F(x-j;\theta_0) + \bar{F}(x;\theta_0)= 1 - \sum_{i=x-j+1}^{x-1} f(i;\theta_0),\] where \(F(x;\theta) = \sum_{i=0}^{x} f(i;\theta)\) and \(\bar{F}(x;\theta) = \sum_{i=x}^{\infty} f(i;\theta) = 1- \sum_{i=0}^{x-1} f(i;\theta)\), with \(F(x;\theta)=0\) when \(x<0\).
Thus, \[p_m(x;\theta_0) =
\left\{ \begin{array}{ll}
F(x-j;\theta_0) + \bar{F}(x;\theta_0) & \mbox{ if }
\tilde{\theta}_{x-j,x} \leq \theta_0 < \tilde{\theta}_{x-j+1,x}, j
=2,...,x+1 \\
1 & \mbox{ if } \tilde{\theta}_{x-2,x} \leq \theta_0
\leq \tilde{\theta}_{x,x+1} \\
F(x;\theta_0) + \bar{F}(x+j+1;\theta_0) & \mbox{ if }
\tilde{\theta}_{x,x+j} < \theta_0 \leq \tilde{\theta}_{x,x+j+1}, j
\geq 1 \\
\end{array}
\right.\] Consider the piecewise continuous part where \(\tilde{\theta}_{x-j,x} \leq \theta_0 <
\tilde{x-j+1,x}\). At \(\theta_0=
\tilde{\theta}_{x-j,x}\), then \(f(x-j;\theta_0)=f(x;\theta_0)\), as \(\theta_0\) increases then \(f(x-j;\theta_0)<f(x;\theta_0)\) and no
other value is added to the \(F(x-j;\theta_0)\) part of the p-value until
\(\theta_0=\tilde{x-j+1,x}\).
The p-value is piecewise continuous part where \(\tilde{\theta}_{x,x+j} < \theta_0 \leq
\tilde{x,x+j+1}\). At \(\theta_0=
\tilde{\theta}_{x,x+j+1}\), then \(f(x;\theta_0)=f(x+j+1;\theta_0)\), as \(\theta_0\) decreases then \(f(x;\theta_0)>f(x+j+1;\theta_0)\) and no
other value is added to the \(\bar{F}(x+j;\theta_0)\) part of the p-value
until \(\theta_0=\tilde{x,x+j}\).
Using derivatives, we can show that during each piecewise continuous
part of the p-value when \(\theta_0>x\), the p-value function is
increasing in \(\theta_0\).
Consider the piecewise continuous part with \(p_m= F(a;\theta_0) + \bar{F}(b;\theta_0)\)
with \(a<b\). First, rewrite \(p_m\) as \[p_m =
\sum_{i=0}^{a} f(i;\theta_0) + 1 - \sum_{j=0}^{b-1}
f(j;\theta_0)\] which simplifies to \[p_m = 1 - \sum_{j=a+1}^{b-1}
f(j;\theta_0).\] Now we take the derivative with respect to \(\theta_0\). Notice that \[ \frac{\partial f(j;\theta)}{\partial \theta} =
\frac{1}{j!} \frac{\partial \theta^j e^{-\theta}}{\partial \theta} =
\frac{1}{j!} \left( j \theta^{j-1} e^{-\theta} - \theta^j e^{-\theta}
\right)= f(j-1;\theta) - f(j;\theta)\] So that \[\frac{ \partial p_m}{\partial \theta_0} = \frac{
\partial \left( 1 - \sum_{j=a+1}^{b-1} f(j;\theta_0) \right) }{\partial
\theta_0}=
- \sum_{j=a+1}^{b-1} \frac{ \partial f(j;\theta_0)}{\partial \theta_0} =
\sum_{j=a+1}^{b-1} \left\{ f(j;\theta_0) - f(j-1;\theta_0) \right\} =
f(b-1;\theta_0) - f(a; \theta_0).\]
When \(a=x\) and \(b=x+j+1\), then when \(\tilde{\theta}_{x,x+j} < \theta_0 \leq \tilde{\theta}_{x,x+j+1}\) and the derivative is \(f(x+j;\theta_0) - f(x;\theta_0).\) Just below the lower end of the interval, \(f(x+j;\theta_0) =f(x;\theta_0)\), then as \(\theta_0\) increases \(f(x+j;\theta_0) >f(x;\theta_0)\) so that the derivative is always positive, and the p-value is increasing.
So the upper confidence limit of the \(100(1-\alpha)\)% minlike CI is the largest \(j\) such that \(p_m(x,\tilde{\theta}_{x,x+j+1})> \alpha\).
For the lower confidence limit we conjecture that the piecewise continuous p-value function is monotonic within each continuous piece.
We compare our results to those of from the function of Reiczigel as referenced in Kaschka and Reiczigel (2020). The function was accessed in May 2021 at the URL referenced in Kaschka and Reiczigel (2020); however, the website was under re-construction in June 2021. I suggest going to https://www2.univet.hu/ and following links to J. Reiczigel’s homepage to find the function (or look at the Rmd file that created this vignette).
x | conf level | lower CL | upper CL | |
---|---|---|---|---|
Reiczigel | 1 | 0.70 | 0.3566749 | 3.309751 |
Fay | 1 | 0.70 | 0.3566749 | 3.309751 |
Reiczigel | 0 | 0.95 | 0.0000000 | 3.764351 |
Fay | 0 | 0.95 | 0.0000000 | 3.764351 |
Reiczigel | 8 | 0.99 | 2.9061062 | 18.807287 |
Fay | 8 | 0.99 | 2.9061062 | 18.807287 |
Reiczigel | 8 | 0.90 | 4.5319366 | 14.512948 |
Fay | 8 | 0.90 | 4.5319366 | 14.512948 |
Reiczigel | 8 | 0.01 | 8.0000000 | 9.000000 |
Fay | 8 | 0.01 | 8.0000000 | 9.000000 |
Blaker’s p-values are (see e.g., Klaschka and Reiczigel, 2020, Section 3.1) \[p_B(x;\theta_0) = \left\{ \begin{array}{ll} F(x; \theta_0) + \bar{F}(x_R; \theta_0) & \mbox{ if $F(x;\theta_0) < \bar{F}(x;\theta_0)$} \\ 1 & \mbox{ if $F(x;\theta_0) = \bar{F}(x;\theta_0)$} \\ F(x_L; \theta_0) + \bar{F}(x; \theta_0) & \mbox{ if $F(x;\theta_0) > \bar{F}(x;\theta_0)$} \\ \end{array} \right. ,\] where \(x_R = min \left\{ i: \bar{F}(i;\theta_0) \leq F(x; \theta_0) \right\}\) and \(x_L = max \left\{ i: {F}(x;\theta_0) \leq \bar{F}(i; \theta_0) \right\}.\)
This p-values are also piecewise continuous, but the jumps are at different places.
Here is a comparison of the three ways of calculating the 95% two-sided p-values for \(x=8\): the central method, the minlike method, and Blaker’s method.
An important computational point is that Blaker’s p-values are not monotonic within each piecewise constant interval. To see this, we focus in on a two pieces of the previous graph.
We conjecture that on the lower side (p-values with \(\theta_0< x\)) the piecewise continuous function is monotonicly increasing within each interval, while on the upper side (p-values with \(\theta_0>x\)) the piecewise continuous function decreases to a inflection point, then increases (see the graph). We can solve for the inflection point by taking the derivative of the p-value function.
For \(a<b\), let \(\hat{\theta}_{a,b}\) be the value \(\theta_0\) such that \[F(a;\theta_0) = \bar{F}(b;\theta_0).\] As
with the ‘minlike’ method, for calculating the confidence
intervals
it is useful to rewrite the p-value function in terms of the piecewise
continuous
functions with jumps at \(\hat{\theta}_{a,b}\) values.
\[p_B(x;\theta_0) = \left\{ \begin{array}{ll} F(x-j; \theta_0) + \bar{F}(x; \theta_0) & \mbox{ if $\hat{\theta}_{x-j,x} \leq \theta_0 < \hat{\theta}_{x-j+1,x}$} \\ 1 & \mbox{ if $\hat{\theta}_{x, x+1} \leq \theta_0 \leq \hat{\theta}_{x-1, x}$} \\ F(x; \theta_0) + \bar{F}(x+j; \theta_0) & \mbox{ if $\hat{\theta}_{x,x+j} < \theta_0 \leq \hat{\theta}_{x,x+j+1}$} \\ \end{array} \right. .\]
Using the results of the previous section, we see that the inflection point on the interval betweeen \(\hat{\theta}_{x,x+j}\) and \(\hat{\theta}_{x,x+j+1}\) is at \(\tilde{\theta}_{x,x+j}\). Using these assumptions, we can solve for the confidence intervals at either the jump points or using uniroot function when necessary.
Klaschka developed the BlakerCI R package, so we can check our results against that.
x | conf level | lower CL | upper CL | |
---|---|---|---|---|
Klaschka | 0 | 0.7000 | 0.0000000 | 1.568120 |
Fay | 0 | 0.7000 | 0.0000000 | 1.568120 |
Klaschka | 1 | 0.9500 | 0.0512933 | 5.525705 |
Fay | 1 | 0.9500 | 0.0512933 | 5.525705 |
Klaschka | 8 | 0.9442 | 3.5501406 | 15.553791 |
Fay | 8 | 0.9442 | 3.5501406 | 15.553791 |
Klaschka | 8 | 0.9440 | 3.5501406 | 15.199944 |
Fay | 8 | 0.9440 | 3.5501406 | 15.199944 |
Klaschka | 8 | 0.9430 | 3.5501406 | 15.118012 |
Fay | 8 | 0.9430 | 3.5501406 | 15.118014 |
Klaschka, J, and Reiczigel, J (2020). On matching confidence intervals and tests for some discrete distributions: methodological and computational aspects. Computational Statistics. DOI: 10.1007/s00180-020-00986-0.