#  Estimating odds and log odds with guaranteed accuracy

Luis Mendo
Universidad Politécnica de Madrid, Spain
<EMAIL_ADDRESS>

###### Abstract

Two sequential estimators are proposed for the odds $p/(1-p)$ and log odds
$\log(p/(1-p))$ respectively, using independent Bernoulli random variables
with parameter $p$ as inputs. The estimators are unbiased, and guarantee that
the variance of the estimation error divided by the true value of the odds, or
the variance of the estimation error of the log odds, are less than a target
value for any $p\in(0,1)$. The estimators are close to optimal in the sense of
Wolfowitz’s bound.

_Keywords:_ Estimation, odds, log odds, mean-square error, efficiency.

_MSC2010:_ 62F10, 62L12.

##  1 Introduction

Consider a probability $p\in(0,1)$, and let $\gamma=p/(1-p)$ and
$\phi=\log\gamma$ respectively denote the _odds_ and _log odds_ associated
with $p$. These are commonly used measures of the probability of an event,
which find application in many fields including clinical trials, epidemiology,
and data analysis. Also, log odds are inherent to logistic regression, which
is a prevalent tool in medical and social sciences (Agresti, 2002; Norton and
Dowd, 2018), and is widely used in machine learning (Bishop, 2006).

For an estimator $\hat{\gamma}$ of $\gamma$, a usual measure of accuracy is
the _mean-square error_ (MSE), $\operatorname{E}[(\hat{\gamma}-\gamma)^{2}]$,
or its square root, known as root-mean-square error (RMSE). Either of these
measures of error is most meaningful in a _relative_ sense. For example, an
RMSE equal to $0.5$ may be adequate if the true odds are $\gamma=9$
(corresponding to $p=0.9$), but unacceptable if $\gamma=0.25$ ($p=0.2$). It is
reasonable to require that the RMSE be proportional to the true (unknown)
value of $\gamma$; or equivalently that the MSE be proportional fo
$\gamma^{2}$. The estimator $\hat{\gamma}$ should ideally guarantee
$\operatorname{E}[(\hat{\gamma}-\gamma)^{2}]/\gamma^{2}$ smaller than a
_prescribed_ value, irrespective of $\gamma$.

When estimating the log odds $\phi$, on the other hand, the error should not
be understood in a relative sense, because the logarithm already has a
normalizing effect, as it transforms ratios into differences. Specifically, a
difference $\delta$ between $\phi$ and an estimate $\hat{\phi}$ thereof
corresponds to a deviation by a factor $\exp\delta$ between $\gamma$ and its
estimate $\hat{\gamma}$ defined as $\exp\hat{\phi}$. That is, relative error
in estimating $\gamma$ corresponds to absolute (non-normalized) error in
estimating $\phi$.

This paper deals with the problem of estimating $\gamma$ and $\phi$ from
independent, identically distributed (i.i.d.) observations of a Bernoulli
distribution with parameter $p$. Sequential estimators of $\gamma$ and $\phi$
are presented that achieve a specified level of accuracy, defined as relative
or absolute MSE respectively. The sequential character of the estimators
cannot be avoided, because a fixed sample size cannot satisfy these conditions
for all $p$. The proposed estimators are based on _inverse binomial sampling_
Haldane, 1945; Mikulski and Smith, 1976; Lehmann and Casella, 1998, chapter 2,
and have a number of desirable features: they are simple to compute; they are
unbiased; and their performance is close to optimal, in the sense that the
relationship between accuracy and average sample size is only slightly worse
than that given by Wolfowitz’s bound (Wolfowitz, 1947).

The problems of estimating the odds and the log odds, as well as related
parameters such as the ratio of two probabilities or the corresponding odds
ratio, have drawn a great deal of attention over the years. However, to the
author’s knowledge, no estimator has been proposed that can guarantee a
specified level of accuracy as discussed above, i.e. in a relative sense for
the odds or in an absolute sense for the log odds, for all values of $p$. A
review of existing results in this area is given next.

Early works by Gart and Zweifel (1967) and Wells and Donner (1987) focus on
analysing the bias and mean-square error of several estimators of the odds
under binomial (i.e. fixed-size) sampling. Siegmund (1982) studies the
asymptotic properties of estimators of the odds and the odds ratio.
Sungboonchoo et al. (2021) propose estimators of the odds with exponentially
decreasing bias, using either binomial or inverse binomial sampling; and
provide asymptotic confidence intervals for estimators of the odds ratio.
Ngamkham et al. (2016) and Ngamkham (2018) consider interval estimation of the
ratio of two probabilities, when each of the two populations is sampled
independently using either a fixed sample size or inverse binomial sampling,
and give asymptotic expressions for confidence intervals with a desired
coverage probability. Similarly, Tian et al. (2008) discuss the estimation of
the ratio of two probabilities using independent inverse binomial sampling for
each population, and review several approaches for obtaining approximate
confidence intervals. Cho (2007) estimates the probability ratio under a loss
function defined as the sum of squared error and cost proportional to the
number of observations, where samples are taken in pairs, one from each
population; and proposes a sequential procedure which is asymptotically
optimal when the cost per sample is small. Cho (2013) presents a sequential
estimator for a probability ratio when the proportion of sample sizes is
specified, and studies its asymptotic properties. Kokaew et al. (2021, 2023)
consider correlation between observations of the two populations, assuming
sample sets of fixed size; and propose estimators of the probability ratio or
its logarithm, for which they derive asymptotic confidence intervals. Agresti
(1999) analyses confidence intervals for the log odds ratio when independent
binomial sampling is applied to each population, and Bandyopadhyay et al.
(2017) address the more general scenario where the two populations are sampled
using different combinations of binomial and inverse binomial sampling.

The following terminology and simple facts will be used throughout the paper.
The two possible values of the Bernoulli distribution will be referred to as
“success” and “failure”, which occur with probabilities $p$ and $1-p$
respectively. In inverse binomial sampling, given $r\in\mathbb{N}$, a sequence
of i.i.d. Bernoulli variables with parameter $p$ is observed until the $r$-th
success is obtained. The required number of samples, $N$, has a _negative
binomial_ distribution:

$\Pr[N=n]=\binom{n-1}{r-1}p^{r}(1-p)^{n-r},\ \ \ n\geq r.$ (1)

(An alternative definition, not used in this paper, considers the negative
binomial distribution as that of $N-r$.) For this distribution Haldane, 1945;
Ross, 2010, chapter 4,

$\displaystyle\operatorname{E}[N]$ $\displaystyle=\frac{r}{p},$ (2)
$\displaystyle\operatorname{E}\left[\frac{1}{N-1}\right]$
$\displaystyle=\frac{p}{r-1},$ (3) $\displaystyle\operatorname{Var}[N]$
$\displaystyle=\frac{r(1-p)}{p^{2}}.$ (4)

##  2 Estimation of odds

The estimator of $\gamma=p/(1-p)$ to be proposed is based on estimating $p$
and $1/(1-p)$ separately by means of inverse binomial sampling. Specifically,
given $r_{1}\in\mathbb{N}$, $r_{1}\geq 3$, a random number $N_{1}$ of samples
are observed until $r_{1}$ _successes_ occur. Then
$\hat{p}_{1}=(r_{1}-1)/(N_{1}-1)$ is an unbiased estimator of $p$ as a
consequence of (3), and its variance satisfies the bound (Pathak and Sathe,
1984)

$\operatorname{Var}[\hat{p}_{1}]\leq\frac{p^{2}(1-p)}{r_{1}-2+2p}.$ (5)

Similarly, for $r_{2}\in\mathbb{N}$, a second, independent set of samples is
taken until $r_{2}$ _failures_ are observed. The number of samples $N_{2}$ is
a negative binomial random variable with parameters $r_{2}$ and $1-p$.
Replacing $r$ and $p$ in (2) and (4) by $r_{2}$ and $1-p$, it follows that
$\hat{d}_{2}=N_{2}/r_{2}$ is an unbiased estimator of $1/(1-p)$, and

$\operatorname{Var}[\hat{d}_{2}]=\frac{p}{r_{2}(1-p)^{2}}.$ (6)

Given $r\in\mathbb{N}$, $r\geq 2$, the estimator of $\gamma$ is defined as

$\hat{\gamma}=\hat{p}_{1}\hat{d}_{2}$ (7)

where $\hat{p}_{1}$ and $\hat{d}_{2}$ are obtained as above with $r_{1}=r+1$
and $r_{2}=r-1$. Let $N=N_{1}+N_{2}$ denote the total number of samples used
by this estimator.

######  Theorem 1.

For $r\in\mathbb{N}$, $r\geq 2$, $p\in(0,1)$, the estimator $\hat{\gamma}$
given by (7) has the following properties:

$\displaystyle\operatorname{E}[\hat{\gamma}]$ $\displaystyle=\gamma$ (8)
$\displaystyle\frac{\operatorname{Var}[\hat{\gamma}]}{\gamma^{2}}$
$\displaystyle\leq\frac{1}{r-1}\left(1-\frac{p(1-p)}{r-1+2p}\right)$ (9)
$\displaystyle<\frac{1}{r-1}$ (10) $\displaystyle\operatorname{E}[N]$
$\displaystyle=\frac{r+1-2p}{p(1-p)}.$ (11)

###### Proof.

The estimators $\hat{p}_{1}$ and $\hat{d}_{2}$ are statistically independent,
and thus

$\operatorname{E}[\hat{\gamma}]=\operatorname{E}[\hat{p}_{1}]\operatorname{E}[\hat{d}_{2}]=\frac{p}{1-p}=\gamma,$
(12)

which proves (8).

Similarly, using (5) and (6) with $r_{1}=r+1$, $r_{2}=r-1$,

$\begin{split}\operatorname{Var}[\hat{\gamma}]&=\operatorname{E}[\hat{p}_{1}^{2}]\operatorname{E}[\hat{d}_{2}^{2}]-\frac{p^{2}}{(1-p)^{2}}\\\
&=\left(\operatorname{Var}[\hat{p}_{1}]+p^{2}\right)\left(\operatorname{Var}[\hat{d}_{2}]+\frac{1}{(1-p)^{2}}\right)-\frac{p^{2}}{(1-p)^{2}}\\\
&\leq\left(\frac{p^{2}(1-p)}{r-1+2p}+p^{2}\right)\left(\frac{p}{(r-1)(1-p)^{2}}+\frac{1}{(1-p)^{2}}\right)-\frac{p^{2}}{(1-p)^{2}}\\\
&=\frac{p^{2}}{(1-p)^{2}}\left(\frac{p(1-p)}{(r-1+2p)(r-1)}+\frac{1-p}{r-1+2p}+\frac{p}{r-1}\right)\\\
&=\frac{(r-1+p+p^{2})p^{2}}{(r-1)(r-1+2p)(1-p)^{2}}\\\
&=\frac{p^{2}}{(r-1)(1-p)^{2}}\left(1-\frac{p(1-p)}{r-1+2p}\right)<\frac{p^{2}}{(r-1)(1-p)^{2}},\end{split}$
(13)

as in (9) and (10).

The average number of samples is obtained as

$\operatorname{E}[N]=\operatorname{E}[N_{1}]+\operatorname{E}[N_{2}]=\frac{r+1}{p}+\frac{r-1}{1-p}=\frac{r+1-2p}{p(1-p)},$
(14)

in agreement with (11). ∎

From Theorem 1 it stems that a relative error as low as desired can be
guaranteed, irrespective of $p$, by choosing $r$ large enough. For example, to
achieve an RMSE smaller than $20\%$ of the true value of $\gamma$ it is
sufficient to take $r=26$, according to (10). Guaranteeing better accuracy
values requires larger $r$, which in turn causes the average sample size to
increase, as expressed by (11).

To assess how good the trade-off between accuracy and average sample size is,
it is meaningful to use Wolfowitz’s bound, which is a generalization of the
Cramér-Rao inequality for sequential estimators Wolfowitz, 1947; Ghosh et al.,
1997, section 4.3. The bound indicates the lowest variance that an estimator
can have for a given average sample size, provided that certain (mild)
regularity conditions are satisfied by the distribution function of the
observations, the sequential procedure and the estimator. When the
observations are i.i.d. Bernoulli random variables with parameter $p$ and the
parameter to be estimated, $\theta$, is a differentiable function of $p$,
Wolfowitz’s bound for an unbiased estimator $\hat{\theta}$ is (DeGroot, 1959,
lemma 2.7)

$\operatorname{Var}[\hat{\theta}]\operatorname{E}[N]\geq\left(\displaystyle\frac{\mathrm{d}\theta}{\mathrm{d}p}\right)^{2}p(1-p).$
(15)

The _efficiency_ of $\hat{\theta}$, which will be denoted as
$\eta(\hat{\theta})$, can then be defined as Wolfowitz’s bound for
$\operatorname{Var}[\hat{\theta}]\operatorname{E}[N]$ divided by the actual
value of that product achieved by $\hat{\theta}$:

$\eta(\hat{\theta})=\frac{({\mathrm{d}\theta}/{\mathrm{d}p})^{2}\,p(1-p)}{\operatorname{Var}[\hat{\theta}]\operatorname{E}[N]}.$
(16)

Taking $\theta=\gamma$, $\hat{\theta}=\hat{\gamma}$, which implies that
$\mathrm{d}\theta/\mathrm{d}p=1/(1-p)^{2}$ and that
$\operatorname{Var}[\hat{\gamma}]$ and $\operatorname{E}[N]$ satisfy (9) and
(11) respectively, (16) becomes

$\displaystyle\eta(\hat{\gamma})$
$\displaystyle=\frac{p}{(1-p)^{3}\operatorname{Var}[\hat{\gamma}]\operatorname{E}[N]}$
(17)
$\displaystyle>\frac{r-1}{r+1-2p}\left(1+\frac{p(1-p)}{r-1+p+p^{2}}\right)$
(18) $\displaystyle>\frac{r-1}{r+1}.$ (19)

Thus, the efficiency of $\hat{\gamma}$ approaches $1$ for large $r$.

Figure 1 compares, for several values of $r$, the lower bound (18) with the
actual efficiency, obtained by Monte Carlo simulation. Specifically, for each
$p$, $10^{8}$ values of the estimate $\hat{\gamma}$ are simulated using
Bernoulli random variables with parameter $p$; the sample variance of
$\hat{\gamma}$ is computed; and this is used instead of the true variance in
(17) to obtain $\eta(\hat{\gamma})$. As can be seen, the bound is quite tight,
especially for $p$ near $0$ or $1$ or for $r$ not too small.

Figure 1: Efficiency of odds estimator $\hat{\gamma}$

##  3 Estimation of log odds

The estimator of $\phi=\log(p/(1-p))$, to be defined precisely later, will be
based on estimating $\log p$ and $\log(1-p)$ separately using inverse binomial
sampling. Thus the estimation of $\log p$ is analysed next.

Given $r\in\mathbb{N}$, $r\geq 2$, and denoting by $N$ the number of samples
required to obtain $r$ successes, an unbiased estimator of $\log p$ can be
obtained from $r$ and $N$ using the following result by DeGroot (1959, theorem
4.1). Let $q=1-p$. Given a function $g(q)$, there exists an unbiased estimator
of $\theta=g(q)$ in inverse binomial sampling if and only if $g(q)$ has a
Taylor expansion around $q=0$ valid for $|q|<1$; and in that case the
estimator $\hat{\theta}(N)$ is given by

$\displaystyle\hat{\theta}(N)$
$\displaystyle=\frac{(r-1)!\,c_{r,N-r}}{(N-1)!},$ (20) $\displaystyle c_{r,k}$
$\displaystyle=\frac{\mathrm{d}^{k}}{\mathrm{d}q^{k}}\left.\frac{g(q)}{(1-q)^{r}}\right\rfloor_{q=0}.$
(21)

Let the $n$-th harmonic number be denoted as $H_{n}=1+1/2+\cdots+1/n$.

######  Proposition 1.

For $g(q)=\log(1-q)=\log p$, (20) and (21) yield

$\hat{\theta}(N)=-\left(\frac{1}{r}+\frac{1}{r+1}+\cdots+\frac{1}{N-1}\right)=-H_{N-1}+H_{r-1}.$
(22)

###### Proof.

Let $x^{(k)}$ denote the rising factorial, $x^{(k)}=x(x+1)\cdots(x+k-1)$,
$x^{(0)}=1$. Direct computation using (21) with $g(q)=\log(1-q)$ shows that
$c_{r,0}=0$, and for $k>0$ this recurrence relation holds:

$c_{r,k}=-r^{(k-1)}+(r+k-1)c_{r,k-1}.$ (23)

From (20) and (23) it stems that $\hat{\theta}(r)=0$, and for $n>r$

$\begin{split}\hat{\theta}(n)&=\frac{(r-1)!\,(-r^{(n-r-1)}+(n-1)c_{r,n-r-1})}{(n-1)!}\\\
&=-\frac{1}{n-1}+\frac{(r-1)!\,c_{r,n-r-1}}{(n-2)!}=-\frac{1}{n-1}+{\hat{\theta}(n-1)},\end{split}$
(24)

which is equivalent to (22). ∎

By Proposition 1, an unbiased estimator $\hat{\lambda}$ of $\log p$ using
inverse binomial sampling is

$\hat{\lambda}=-H_{N-1}+H_{r-1},$ (25)

where $r$ is the number of successes and $N\geq r$ is the number of
observations. The variance of this estimator is bounded as shown next.

######  Theorem 2.

For $r\in\mathbb{N}$, $r\geq 2$, $p\in(0,1)$, the variance of $\hat{\lambda}$
satisfies

$\displaystyle\operatorname{Var}[\hat{\lambda}]$
$\displaystyle<\frac{1}{r-1+p}\left(\left(1+\frac{p}{2r-1}-\frac{1}{4r(2r-1)}\right)(1-p)+\frac{p}{4(r-1)}\right)$
(26) $\displaystyle<\frac{1}{r-1}.$ (27)

The proof of this theorem makes use of the following lemma.

######  Lemma 1.

For $r,\,n\in\mathbb{N}$, $n>r$,

$\log\frac{n-1/2+\alpha}{r-1/2+\alpha}<H_{n-1}-H_{r-1}<\log\frac{n-1/2}{r-1/2}$
(28)

with

$\displaystyle\alpha$ $\displaystyle=\frac{-r+\sqrt{r^{2}+1}}{2},$ (29)
$\displaystyle 0<\alpha$ $\displaystyle<\frac{1}{4r}.$ (30)

###### Proof.

Consider $r,\,n\in\mathbb{N}$, $n\geq r$. The shaded area in Figure 2(a),
where the heights of the rectangles are $1/r,\ldots,1/(n-1)$, equals
$H_{n-1}-H_{r-1}$. Replacing the top side of each rectangle by the straight
line that passes through its middle point and is tangent to the curve
$y=1/(x+1/2)$, a trapezium is obtained with the same area as the rectangle.
The figure shows, with dashed line, the trapezium associated with one of the
rectangles. It follows from the convexity of the curve that

$H_{n-1}-H_{r-1}<\int_{r-1}^{n-1}\frac{\mathrm{d}x}{x+1/2}=\log\frac{n-1/2}{r-1/2},$
(31)

which establishes the second inequality in (28).

(a) Upper bound

(b) Lower bound

Figure 2: Proof of bounds in (28)

For the first inequality in (28), consider the shaded rectangle in Figure
2(b), which extends from $x=k-1$ to $x=k$ and has height $1/k$. Let
$\alpha_{k}$ be chosen so that the intercept points of the curve
$y=1/(x+1/2+\alpha_{k})$ with the vertical lines $x=k-1$ and $x=k$ are
vertically equispaced from the horizontal line $y=1/k$, i.e. have heights
$1/k+\beta_{k}$ and $1/k-\beta_{k}$ for some $\beta_{k}>0$. That is,
$\alpha_{k}$ and $\beta_{k}$ satisfy the equation system

$\displaystyle\frac{1}{k-1/2+\alpha_{k}}$
$\displaystyle=\frac{1}{k}+\beta_{k}$ (32)
$\displaystyle\frac{1}{k+1/2+\alpha_{k}}$
$\displaystyle=\frac{1}{k}-\beta_{k},$ (33)

which yields

$\alpha_{k}=\frac{-k+\sqrt{k^{2}+1}}{2}.$ (34)

By construction, the area of the shaded rectangle, which is $1/k$, coincides
with that of the trapezium shown with dashed line in the figure, whose upper
side is the segment that joins the two intercept points of the curve. This
segment is a chord of the curve, which is convex, and therefore the area under
the curve between $x=k-1$ and $x=k$ is smaller than the shaded area:

$\int_{k-1}^{k}\frac{\mathrm{d}x}{x+1/2+\alpha_{k}}<\frac{1}{k}.$ (35)

The terms $\alpha_{r},\ldots,\alpha_{n-1}$ as given by (34) form a decreasing
sequence. Thus defining $\alpha=\alpha_{r}$, as in (29), ensures that
$\alpha>\alpha_{k}$ for $k>r$, and therefore

$\begin{split}\int_{r-1}^{n-1}\frac{\mathrm{d}x}{x+1/2+\alpha}&=\log\frac{n-1/2+\alpha}{r-1/2+\alpha}\\\
&<\sum_{k=r}^{n-1}\int_{k-1}^{k}\frac{\mathrm{d}x}{x+1/2+\alpha_{k}}<\sum_{k=r}^{n-1}\frac{1}{k}=H_{n-1}-H_{r-1},\end{split}$
(36)

which proves the first inequality in (28).

Lastly, rewriting (29) as

$\alpha=\frac{-r+r\sqrt{1+1/r^{2}}}{2}$ (37)

and using the inequality $1<\sqrt{1+x}<1+x/2$, $x\neq 0$ yields (30). ∎

###### Proof of Theorem 2.

Consider $r\in\mathbb{N}$, $r\geq 2$ and $p\in(0,1)$, and let $\alpha$ be
defined as in Lemma 1. Noting that $\log((n-1/2-\zeta)/(r-1/2-\zeta))$ is a
continuous function of $\zeta\in[0,\alpha]$, it stems from the referred lemma
that there is $\zeta_{n}\in(0,\alpha)$ such that

$H_{n-1}-H_{r-1}=\log\frac{n-1/2+\zeta_{n}}{r-1/2+\zeta_{n}}.$ (38)

Therefore, $\operatorname{Var}[\hat{\lambda}]$ can be written as

$\operatorname{Var}[\hat{\lambda}]=\operatorname{E}[(-H_{N-1}+H_{r-1}-\log
p)^{2}]=\operatorname{E}\left[\log^{2}\frac{r-1/2+\zeta_{N}}{(N-1/2+\zeta_{N})p}\right],$
(39)

where $\zeta_{N}$ is a random variable that takes values in $(0,\alpha)$. For
convenience, let

$\alpha^{\prime}=\frac{1}{4r}.$ (40)

Making use of the inequality

$\log^{2}x<\frac{1}{x}+x-2,\ \ \ x>0,$ (41)

which results from Topsøe (2007, equation (14)), and taking into account that
$0<\zeta_{N}<\alpha<\alpha^{\prime}$, (39) becomes

$\displaystyle\operatorname{Var}[\hat{\lambda}]$
$\displaystyle<\operatorname{E}\left[\frac{(N-1/2+\zeta_{N})p}{r-1/2+\zeta_{N}}\right]+\operatorname{E}\left[\frac{r-1/2+\zeta_{N}}{(N-1/2+\zeta_{N})p}\right]-2$
(42)
$\displaystyle<\operatorname{E}\left[\frac{(N-1/2)p}{r-1/2}\right]+\operatorname{E}\left[\frac{r-1/2+\alpha^{\prime}}{(N-1/2+\alpha^{\prime})p}\right]-2.$
(43)

Computing the first summand in (43) is straightforward, in view of (2). For
the second, an upper bound is obtained next. It can be easily checked that,
for $n>1$ and $|\alpha^{\prime}|<1/2$,

$\frac{1}{n-1/2+\alpha^{\prime}}<\frac{1/2-\alpha^{\prime}}{n-1}+\frac{1/2+\alpha^{\prime}}{n}.$
(44)

Using (3) and the inequality (Mendo and Hernando, 2006, section II)

$\operatorname{E}\left[\frac{1}{N}\right]<\frac{p}{r-1+p},$ (45)

it follows from (44) that

$\operatorname{E}\left[\frac{1}{(N-1/2+\alpha^{\prime})p}\right]<\frac{1/2-\alpha^{\prime}}{r-1}+\frac{1/2+\alpha^{\prime}}{r-1+p}=\frac{r-1+(1/2-\alpha^{\prime})p}{(r-1)(r-1+p)}.$
(46)

Using (2) and (46) into (43),

$\begin{split}\operatorname{Var}[\hat{\lambda}]&<\frac{r-p/2}{r-1/2}+\frac{(r-1/2+\alpha^{\prime})(r-1+(1/2-\alpha^{\prime})p)}{(r-1)(r-1+p)}-2\\\
&=\frac{1-p}{2r-1}+\frac{(r-1)(1/2+\alpha^{\prime})(1-p)+(1/4-\alpha^{\prime
2})p}{(r-1)(r-1+p)}\\\
&=\left(\frac{1/2}{r-1/2}+\frac{1/2+\alpha^{\prime}}{r-1+p}\right)(1-p)+\frac{(1/4-\alpha^{\prime
2})p}{(r-1)(r-1+p)}.\end{split}$ (47)

Substituting (40),

$\begin{split}\operatorname{Var}[\hat{\lambda}]&<\left(\frac{1/2}{r-1/2}+\frac{1/2+1/(4r)}{r-1+p}\right)(1-p)+\frac{(1/4-1/(16r^{2}))p}{(r-1)(r-1+p)}\\\
&=\left(\frac{1}{r-1+p}+\frac{p/2-1/(8r)}{(r-1+p)(r-1/2)}\right)(1-p)+\frac{(1/4-1/(16r^{2}))p}{(r-1)(r-1+p)}\\\
&<\frac{1}{r-1+p}\left(\left(1+\frac{p-1/(4r)}{2r-1}\right)(1-p)+\frac{p}{4(r-1)}\right),\end{split}$
(48)

which gives (26).

As a result of the above,

$\displaystyle\operatorname{Var}[\hat{\lambda}]$
$\displaystyle<\frac{1}{r-1}\left(\left(1+\frac{p}{2r-1}\right)(1-p)+\frac{p}{4(r-1)}\right)$
(49)
$\displaystyle<\frac{1}{r-1}\,\frac{2r^{2}(1-p)-r(3-9p/2+p^{2})+1-9p/4+p^{2}}{(2r-1)(r-1)}.$
(50)

The second fraction in (50) is easily seen to be less than $1$, and thus (27)
holds. ∎

Based on Theorem 2, the estimator of $\phi=\log(p/(1-p))$ is defined as
follows. Given $r\in\mathbb{N}$, $r\geq 2$, inverse binomial sampling with $r$
successes is applied to a sequence of i.i.d. Bernoulli random variables
$X_{1i}$, $i=1,2,\ldots$ with parameter $p$. This requires a random number of
samples $N_{1}$, from which an estimate $\hat{\lambda}_{1}$ of $\log p$ is
obtained as in (25). An analogous procedure, with the same $r$, is then
applied to a sequence of “complementary” Bernoulli random variables
$1-X_{2j}$, $j=1,2,\ldots$, where the variables $X_{2j}$ are again i.i.d. with
parameter $p$, and independent of the previously used $X_{1i}$. This consumes
a random number of samples $N_{2}$, and provides an estimate
$\hat{\lambda}_{2}$ of $\log(1-p)$. The estimate $\hat{\phi}$ of $\phi$ is
then obtained as

$\hat{\phi}=\hat{\lambda}_{1}-\hat{\lambda}_{2}=-H_{N_{1}-1}+H_{N_{2}-1}.$
(51)

Let $N=N_{1}+N_{2}$ be the total number of samples used by the estimator
$\hat{\phi}$.

######  Theorem 3.

For $r\in\mathbb{N}$, $r\geq 2$, $p\in(0,1)$, the estimator $\hat{\phi}$
defined by (51) has the following properties:

$\displaystyle\operatorname{E}[\hat{\phi}]$ $\displaystyle=\log\frac{p}{1-p}$
(52) $\displaystyle\operatorname{Var}[\hat{\phi}]$
$\displaystyle<\frac{r^{2}-r/4-1/4}{(r-1+p)(r-p)(r-1/2)}-\frac{p(1-p)}{(r-1/2)^{2}}\left(1-\frac{1}{2r-3}\right)$
(53) $\displaystyle<\frac{1}{r-5/4}\left(1-\frac{7}{(4r-1)^{2}}\right)$ (54)
$\displaystyle\operatorname{E}[N]$ $\displaystyle=\frac{r}{p(1-p)}.$ (55)

###### Proof.

The equality (52) follows from (51) and the fact that $\hat{\lambda}_{1}$ and
$\hat{\lambda}_{2}$ are unbiased.

The bound (26) from Theorem 2 can be rewritten as

$\operatorname{Var}[\hat{\lambda}]<\left(\frac{1}{2}+\frac{1}{4r}\right)\frac{1-p}{r-1+p}+\frac{1-p}{2r-1}+\frac{p}{4(r-1+p)(r-1)}.$
(56)

Since $\hat{\lambda}_{1}$ and $\hat{\lambda}_{2}$ are statistically
independent,
$\operatorname{Var}[\hat{\phi}]=\operatorname{Var}[\hat{\lambda}_{1}]+\operatorname{Var}[\hat{\lambda}_{2}]$.
Applying (56) to bound $\operatorname{Var}[\hat{\lambda}_{1}]$ and
$\operatorname{Var}[\hat{\lambda}_{2}]$, with $p$ replaced by $1-p$ for the
latter, gives

$\begin{split}\operatorname{Var}[\hat{\phi}]&<\left(\frac{1}{2}+\frac{1}{4r}\right)\left(\frac{1-p}{r-1+p}+\frac{p}{r-p}\right)+\frac{1}{2r-1}\\\
&\quad+\frac{1}{4(r-1)}\left(\frac{p}{r-1+p}+\frac{1-p}{r-p}\right)\\\
&=\frac{r-2p(1-p)}{2(r-1+p)(r-p)}+\frac{r-2p(1-p)}{4r(r-1+p)(r-p)}+\frac{1}{2r-1}\\\
&+\frac{1}{4(r-1)}\,\frac{r-1+2p(1-p)}{(r-1+p)(r-p)}.\end{split}$ (57)

Rearranging,

$\begin{split}\operatorname{Var}[\hat{\phi}]&<\frac{1}{2r-1}+\frac{r+1}{2(r-1+p)(r-p)}-\frac{p(1-p)}{(r-1+p)(r-p)}\left(1-\frac{1}{2r(r-1)}\right)\\\
&=\frac{r^{2}-r/4-1/4}{(r-1+p)(r-p)(r-1/2)}\\\
&\quad-\frac{p(1-p)}{(r-1+p)(r-p)}\left(1-\frac{1/2}{r-1/2}-\frac{1/2}{r(r-1)}\right),\end{split}$
(58)

which combined with the inequalities $(r-1+p)(r-p)\leq(r-1/2)^{2}$ and

$\frac{1}{r-1/2}+\frac{1}{r(r-1)}=\frac{r^{2}-1/2}{r(r-1/2)(r-1)}<\frac{r^{2}-1/4}{r(r-1/2)(r-1)}<\frac{1}{r-3/2}$
(59)

produces (53).

Differentiating with respect to $p$, the right-hand side of (53) is seen to be
decreasing for $p<1/2$ and increasing for $p>1/2$. Therefore its supremum
occurs, by symmetry, when $p$ tends to either $0$ or $1$. This gives

$\begin{split}\operatorname{Var}[\hat{\phi}]&<\frac{r^{2}-r/4-1/4}{r(r-1/2)(r-1)}=\frac{r^{3}-3r^{2}/2+r/16+5/16}{r(r-1/2)(r-1)(r-5/4)}\\\
&=\frac{1}{r-5/4}\left(1-\frac{7r-5}{16r(r-1/2)(r-1)}\right)\\\
&<\frac{1}{r-5/4}\left(1-\frac{7}{16r(r-1/2)}\right)<\frac{1}{r-5/4}\left(1-\frac{7}{16(r-1/4)^{2}}\right)\\\
&=\frac{1}{r-5/4}\left(1-\frac{7}{(4r-1)^{2}}\right),\end{split}$ (60)

in accordance with (54).

Lastly, $E[N]=E[N_{1}]+E[N_{2}]=r/p+r/(1-p)=r/(p(1-p))$, which proves (55). ∎

By Theorem 3, $\hat{\phi}$ is an unbiased estimator of $\log(p/(1-p))$.
Furthermore, choosing $r$ large enough guarantees, according to (54), an error
variance as low as desired irrespective of $p$; at the cost of an increase in
average sample size proportional to $r$, as given by (55). Thus, for example,
to ensure an RMSE smaller than $0.2$ it suffices to use $r=27$ (which actually
guarantees that the RMSE is less than $0.19701$).

As in Section 2, the efficiency of $\hat{\phi}$ can be used as a measure of
quality of this estimator. Particularizing (16) for $\theta=\phi$ and
$\hat{\theta}=\hat{\phi}$, which gives
$\mathrm{d}\theta/\mathrm{d}p=1/(p(1-p))$, and making use of (53)–(55),

$\displaystyle\eta(\hat{\phi})$
$\displaystyle=\frac{1}{p(1-p)\operatorname{Var}[\hat{\gamma}]\operatorname{E}[N]}$
(61)
$\displaystyle>\frac{1}{\displaystyle\frac{r^{3}-r^{2}/4-r/4}{(r-1+p)(r-p)(r-1/2)}-\frac{rp(1-p)}{(r-1/2)^{2}}\left(1-\frac{1}{2r-3}\right)}$
(62)
$\displaystyle>\frac{r-5/4}{r}\,\frac{(4r-1)^{2}}{(4r-1)^{2}-7}>\frac{r-5/4}{r}.$
(63)

This shows that, irrespective of $p$, the efficiency of $\hat{\phi}$ is close
to $1$ for moderate or large values of $r$, and tends to $1$ as $r$ grows.

The lower bound (62) is compared in Figure 3 with efficiency values obtained
by Monte Carlo simulation, using the same procedure as in Section 2. As can be
seen, the bound is less tight in this case, compared with the bound (18) for
$\eta(\hat{\gamma})$. This could be expected, because the expression of
$\hat{\lambda}$, on which $\hat{\phi}$ is based, is more difficult to deal
with than that of $\hat{\gamma}$ (in the proof of Theorem 2,
$\operatorname{Var}[\hat{\lambda}]$ is expressed as a series involving the
square of a certain logarithm, which is then bounded by a simpler function).
Still, the uniform bound (63) is greater than that for $\eta(\hat{\gamma})$,
given by (19); and, as was the case for that estimator, it becomes tighter as
$r$ increases.

Figure 3: Efficiency of log odds estimator $\hat{\phi}$

##  4 Discussion and further work

Two sequential estimators have been proposed to estimate the odds $p/(1-p)$
and the log odds $\log(1/(1-p))$, respectively, from Bernoulli observations.
The estimators are unbiased, and guarantee that the relative variance (for
odds) or the variance (for log odds) are below a desired value, irrespective
of $p$. The guaranteed accuracy level is controlled by the parameter $r$.
Larger values of $r$ guarantee better accuracy levels, at the cost of
increased average sample size. The two estimators are close to optimal in the
sense of Wolfowitz’s bound, with efficiency values greater than $(r-1)/(r+1)$
and $(r-5/4)/r$ respectively.

The estimators are based on applying inverse binomial sampling separately with
parameters $p$ and $1-p$, where the latter is obtained by swapping the role of
successes and failures in the Bernoulli samples. This approach might appear
wasteful, because information about both $p$ and $1-p$ could in principle be
inferred using only one of the two sample sets. However, since the two
probabilities are complementary, it is never the case that the two separate
procedures require a large number of samples at the same time; and indeed
performance turns out to be close to optimal, as stated.

On the other hand, using a single set of samples to infer $p/(1-p)$ or its
logarithm would pose some problems. As shown by DeGroot (1959, theorem 4.1),
there is no unbiased estimator of these parameters in inverse binomial
sampling. For $(1-p)/p$ an unbiased estimator does exist; but then it is not
possible to guarantee a relative MSE, because as $p$ grows the sample size
decreases, when it should increase so as to reduce the MSE. By similar
arguments, an unbiased estimator of $\log((1-p)/p)$ exists but is not useful
either. For fixed-size sampling there are no unbiased estimators of any of
these quantities (DeGroot, 1959, section 4); and in any case, with a fixed
sample size it is not possible to guarantee a desired relative MSE for odds,
or a desired MSE for log odds. Moreover, regardless of the rule used for
determining sample size, using a single sample set can make the error analysis
more difficult. Namely, if the estimator is expressed as the product or sum of
two estimators (as has been done in Sections 2 and 3), those two estimators
are not statistically independent, because they are based on the same
observations.

The approach used in this paper could in principle be extended to estimating
the probability ratio (or relative risk) $p_{1}/p_{2}$ or the odds ratio
$p_{1}(1-p_{2})/(p_{2}(1-p_{1}))$, as well as their logarithms, for two
populations with parameters $p_{1}$ and $p_{2}$. For instance, to estimate
$p_{1}/p_{2}$ with guaranteed relative MSE, two inverse binomial sampling
procedures would be used for estimating $p_{1}$ and $1/p_{2}$ respectively;
and error bounds analogous to those derived in the paper are easy to obtain.
However, an important limitation of using this method with two populations is
that it does not provide control on the proportion of samples taken from each
population. Although it is possible to choose different pairs of $r_{1}$,
$r_{2}$ that guarantee approximately the same estimation accuracy, the choice
of $r_{1}$, $r_{2}$ does not determine the distribution of $N_{1}/N_{2}$, or
the ratio $\operatorname{E}[N_{1}]/\operatorname{E}[N_{2}]$, as these depend
also on $p_{1}$, $p_{2}$. Thus, each value of the unknown ratio $p_{1}/p_{2}$
would require a different choice of $r_{1}$, $r_{2}$ in order to attain a
given $\operatorname{E}[N_{1}]/\operatorname{E}[N_{2}]$ and a desired
accuracy.

Having control on $N_{1}/N_{2}$, or at least on
$\operatorname{E}[N_{1}]/\operatorname{E}[N_{2}]$, when sampling from two
populations is often important in applications. For example, if observations
are taken in pairs, one sample from each population, then necessarily
$N_{1}=N_{2}$. Further work is necessary to devise a method that can
incorporate this type of restriction while guaranteeing a given accuracy.

## References

  * Agresti (1999) Alan Agresti.  On logit confidence intervals for the odds ratio with small samples.  _Biometrics_ , 5(2):597–602, June 1999.
  * Agresti (2002) Alan Agresti.  _Categorical Data Analysis_.  John Wiley and Sons, second edition, 2002.
  * Bandyopadhyay et al. (2017) Uttam Bandyopadhyay, Suman Sarkar, and Atanu Biswas.  Fixed-width confidence interval of log odds ratio for joint binomial and inverse binomial sampling.  _Sequential Analysis_ , 36(3):345–354, 2017.
  * Bishop (2006) Christopher M. Bishop.  _Pattern Recognition and Machine Learning_.  Springer, second edition, 2006.
  * Cho (2007) Hokwon Cho.  Confidence intervals for the risk ratio under inverse sampling.  _Journal of Statistical Planning and Inference_ , 137(7):2336–2356, July 2007.  doi: https://doi.org/10.1016/j.jspi.2006.08.005.
  * Cho (2013) Hokwon Cho.  Approximate confidence limits for the ratio of two binomial variates with unequal sample sizes.  _Communications for Statistical Applications and Methods_ , 20(5):347–356, 2013.  doi: https://dx.doi.org/10.5351/CSAM.2013.20.5.347.
  * DeGroot (1959) Morris H. DeGroot.  Unbiased sequential estimation for binomial populations.  _Annals of Mathematical Statistics_ , 30(1):80–101, March 1959.
  * Gart and Zweifel (1967) John J. Gart and James R. Zweifel.  On the bias of various estimators of the logit and its variance with application to quantal bioassay.  _Biometrika_ , 54(1/2):181–187, 1967.
  * Ghosh et al. (1997) Malay Ghosh, Nitis Mukhopadhyay, and Pranab K. Sen.  _Sequential Estimation_.  John Wiley and Sons, 1997.
  * Haldane (1945) J. B. S. Haldane.  On a method of estimating frequencies.  _Biometrika_ , 33(3):222–225, 1945.
  * Kokaew et al. (2021) Angkana Kokaew, Jeeraporn Thaithanan, Winai Bodhisuwan, and Andrei Volodin.  Confidence estimation of a ratio of binomial proportions for dependent populations.  _Lobachevskii Journal of Mathematics_ , 42(2):394–403, 2021.
  * Kokaew et al. (2023) Angkana Kokaew, Winai Bodhisuwan, Su-Fen Yang, and Andrei Volodin.  Logarithmic confidence estimation of a ratio of binomial proportions for dependent populations.  _Journal of Applied Statistics_ , 50(8):1750–1771, 2023.
  * Lehmann and Casella (1998) E. L. Lehmann and George Casella.  _Theory of Point Estimation_.  Springer, second edition, 1998.
  * Mendo and Hernando (2006) Luis Mendo and José María Hernando.  A simple sequential stopping rule for Monte Carlo simulation.  _IEEE Transactions on Communications_ , 54(2):231–241, February 2006.
  * Mikulski and Smith (1976) Piotr W. Mikulski and Paul J. Smith.  A variance bound for unbiased estimation in inverse sampling.  _Biometrika_ , 63(1):216–217, 1976.
  * Ngamkham (2018) Thuntida Ngamkham.  _Confidence Interval Estimation for the Ratio of Binomial Proportions and Random Numbers Generation for some Statistical Models_.  PhD thesis, University of Regina, March 2018.
  * Ngamkham et al. (2016) Thuntida Ngamkham, Andrei Volodin, and Igor N. Volodin.  Confidence intervals for a ratio of binomial proportions based on direct and inverse sampling schemes.  _Lobachevskii Journal of Mathematics_ , 37(4):466–495, 2016.
  * Norton and Dowd (2018) Edward C. Norton and Bryan E. Dowd.  Log odds and the interpretation of logit models.  _Health Services Research_ , 53(2):859–878, April 2018.
  * Pathak and Sathe (1984) P. K. Pathak and Y. S. Sathe.  A new variance formula for unbiased estimation in inverse sampling.  _Sankhyā: The Indian Journal of Statistics, Series B_ , 46(3):301–305, December 1984.
  * Ross (2010) Sheldon Ross.  _A First Course in Probability_.  Pearson Prentice Hall, eighth edition, 2010.
  * Siegmund (1982) David Siegmund.  A sequential confidence interval for the odds ratio.  _Probability and Mathematical Statistics_ , 2(2):149–156, 1982.  doi: https://dx.doi.org/10.5351/CSAM.2013.20.5.347.
  * Sungboonchoo et al. (2021) Chanakan Sungboonchoo, Thuntida Ngamkham, Wararit Panichkitkosolkul, and Andrei Volodin.  Confidence estimation of the cross-product ratio of binomial proportions under different sampling schemes.  _Lobachevskii Journal of Mathematics_ , 42(2):435–450, 2021.
  * Tian et al. (2008) M. Tian, M. L. Tang, H. K. T. Ng, and P. S. Chan.  Confidence intervals for the risk ratio under inverse sampling.  _Statistics in Medicine_ , 27(17):3301–3324, 2008\.  doi: https://doi.org/10.1002/sim.3158.
  * Topsøe (2007) Flemming Topsøe.  Some bounds for the logarithmic function.  In Yeol Je Cho, Jong Kyu Kim, and Sever S. Dragomir, editors, _Inequality Theory and Applications_ , volume 4, pages 137–151. Nova Science Publishers, New York, United States, 2007.  ISBN 13 978-1-59454-874-1.
  * Wells and Donner (1987) G.A. Wells and A. Donner.  Development of formulas for the bias and mean square error of the logit estimator.  In Ian B. MacNeill, Gary J. Umphrey, Allan Donner, and V. Krishna Jandhyala, editors, _Biostatistics: Advances in Statistical Sciences Festschrift in Honor of Professor V.M. Joshi’s 70th Birthday Volume V_ , volume 38, pages 103–124. Springer Netherlands, Dordrecht, 1987.  ISBN 978-94-010-8626-4.  doi: https://doi.org/10.1007/978-94-009-4794-8˙6.
  * Wolfowitz (1947) J. Wolfowitz.  The efficiency of sequential estimates and Wald’s equation for sequential processes.  _Annals of Mathematical Statistics_ , 18(2):215–230, 1947.