First of all I'd like to precise that I'm not an expert of the subject.

Suppose to have two random variables $X$ and $Y$ that are binomial, respectively $X\sim B(n_1,p)$ and $Y\sim B(n_2,p),$ note here that $p$ is the same. I know that $Z=X+Y \sim B(n_1+n_2,p).$

Let $\{x_1,\ldots,x_k\}$ be a sample for $X$ and $\{y_1,\ldots,y_k\}$ be a sample for $Y$, is there a standard method for estimating $n=n_1+n_2$ and $p$?

This is what we have done:

  1. take the "new sample" for $Z$ given by $\{x_1+y_1,\ldots, x_k+y_k\}$,
  2. using the Likelihood Estimator, we obtain estimations for $n$ and $p$,
  3. with the Fisher information, we try to understand the errors over $n$ and $p$.

The method seems to work, but we still have some doubts. Let $S_k$ the group of permutation over $k$ elements. For every $\sigma\in S_k$ we can consider the "sample" given by $\{x_1+y_{\sigma(1)},\dots, x_k+y_{\sigma(k)}\}.$ Applying the Likelihood Estimator for each one of the "new samples" (there are $k!$ different sums) we obtain different estimation $(n_\sigma,p_\sigma)$ for $n$ and $p$ .

What is the meaning of this? How does the new values $n_\sigma, p_\sigma$ are correlated? It can be used for calculating the error for $n$?

Some Comments: The question was previously posted here, but a user suggests me to use stats/crossvalidated SE.

In the example that I have in mind $n$ is the number of birds in a given region and $p$ the visibility probability. I need to aggregate regions with similar $p$, otherwise the data are too small. In particular I need, if possible, an estimation only for $n$, where $p$ a priori is unknown

An example In order to be clear and in view of the answer of kjetil b halvorsen, I will try to put here a practical example. Suppose that we have only one region divided in two zones with probability equal to a fixed $p$ and our data are as follows:

Zone 1   Zone 2
  a1      b1
  a2      b2
  a3      b3
  a4      b4
  a5      b5
  a6      b6

We then can consider this:

Zone 1+2

Then we can use the loglikelihood method in order to estimate $N_1+N_2$ and also $p$ where $N_i$ is the parameter for the binomial of the observed variables in Zone $i$. Is it right?

Now, I know that the likelihood method isn't stable (for me stable means only good). Can we use the Fisher information? If yes, which kind of information can we have?

Finally let be $\sigma$ and $\tau$ two permutations over $6$ elements (there are $(6!)^2$ different couples) than we can consider the new data given by

Zone 1+2
c1=a$\sigma(1)$+ b$\tau(1)$
c2=a$\sigma(2)$+ b$\tau(2)$
c3=a$\sigma(3)$+ b$\tau(3)$
c4=a$\sigma(4)$+ b$\tau(4)$
c5=a$\sigma(5)$+ b$\tau(5)$
c6=a$\sigma(6)$+ b$\tau(6)$

Redoing the likelihood method, with this new variables, we obtain different estimations for $N_1+N_2$.

So the question is: does the set of estimation give me some information about the errors?

  • $\begingroup$ You should give some more detail. How many regions do you want to use (assume has the same visibility parameter $p$? (the Q assumes two))? How many samples do you have for each region? (I assume those are counts over disjoint time intervals of equal length?) Since your interest parameter is $n_1$ and $n_2$ the approach of summing them cannot be correct! will come back when you have answered ... $\endgroup$ Commented Nov 10, 2014 at 15:30
  • $\begingroup$ @kjetilbhalvorsen: The number of regions in my Q is $2$ and I have 6 samples for each region. The counts are over intervals of equal length (but $n_i$ is assumed to be constant). I'm interested in $n_1+n_2$ and not in the single $n_i$. I hope this can help. $\endgroup$ Commented Nov 10, 2014 at 15:40
  • $\begingroup$ Another question: How large (approximately) are the counts? In the tens? hundreds? thousands? Have you any idea about probable values of $p$? Very low? (If the counts are larger, then maybe we can try a poisson approximation? $\endgroup$ Commented Nov 10, 2014 at 15:47
  • $\begingroup$ I tried to aggregate counts because the numbers were low. In reality I have 5 different zone (a zone correspond to regions with the same probability.) In each zone I have something like $15$ regions. The various counts vary from $0$ to $25$ (the $n_i$ seems to be not related and can vary from $0$ to something more than $25$). But the aggregate data can vary from the permutation that we choose. (For example in a zone of $10$ regions with $6$ counts for each region there are $(6!)^{10}$ different aggregations). Finally $p$ seem to be close to $.7$. $\endgroup$ Commented Nov 10, 2014 at 15:56
  • 2
    $\begingroup$ You will also likely benefit a lot from having informative priors, as there is an identification problem between n and p. $\endgroup$
    – Arthur B.
    Commented Nov 13, 2014 at 20:03

2 Answers 2


I will try an answer, even if I am not completely clear about the situation. Formulas will have to be adapted! The problem of estimation of $N$ in the binomial distribution is old, and there are multiple relevant papers. I will give some references at the end.

Let there be $R$ regions (in OP example $R=2$), with $T$ samples (from disjoint time intervals of equal length) from each region. The observed variables is $x_{it}$ which are independent binomial random variables, each with the distribution $\text{Bin}(N_i,p)$ both unknown. The log-likelihood function becomes $$ \ell ( N_i , p ) = \sum \ln \binom{N_i}{x_{it}} + \ln p \cdot \sum x_{it} + \ln (1-p) \cdot \sum ( N_i - x_{it} ) $$ Note that, in the usual problem when $N_i$ is known so that only $p$ is unknown, then the sum (or the mean) of the binomial counts $x_{it}$ is a sufficient summary, so the analysis can be done in terms of the binomial distribution of the sum. In our problem, however, because of the first term in the log-likelihood function, such is not the case, and the log likelihood depends on each of the counts individually! So what you propose, to reduce to the sum of the counts (over $i$), SHOULD NOT BE DONE, as that will lose information (how much, I don't know, but that can be investigated ...). Let us try to understand this a little better. First, we see below that $\max_t(x_{it})$ is a consistent estimator of $N_i$, but this consistent estimator is not a function of the summed counts. That is one clear indication that summation looses information! Note also that the mean is an unbiased estimator of its expectation which is $N_i p$, but doesn't seem to hold information about $N_i$ and $p$ individually, when nothing is known about the other parameter. That indicates that if there is useful information about $N_i$ in the likelihood function, that must be contained in the spread of the values $x_{i1}\dots, x_{iT}$, again indicating that summation is bad. The Olkin et al paper referenced below shows indeed that the method-of-moments estimator in many cases is better than maximum likelihood! and that uses the empirical variance of the $x_{i1}\dots, x_{iT}$, so couldn't be calculated from the summed data.

This problem is known to be unstable. Let us try to understand why. In the usual problem, estimating $p$ when $N_i$ in known, the estimation can be done from some gross characteristic of the data, the mean. When trying to estimate both $N_i$ and $p$, we use much finer properties of the log-likelihood function (thus of the data). To see why, remember that we can obtain the Poisson distribution as a limit of the binomial when $p$ goes to zero and $N$ grows without bounds, with a constant positive product. So, if $p$ is small and $N$ large, the binomial distribution will be quite close to that limit. Take two cases: (A) $N=100, p=0.01$, (B) $N=20, p=0.05$. Draw histograms for the two (binomial) distributions:

zapsmall(cbind(0:20, pA, pB))
               pA       pB
 [1,]  0 0.366032 0.358486
 [2,]  1 0.369730 0.377354
 [3,]  2 0.184865 0.188677
 [4,]  3 0.060999 0.059582
 [5,]  4 0.014942 0.013328
 [6,]  5 0.002898 0.002245
 [7,]  6 0.000463 0.000295
 [8,]  7 0.000063 0.000031
 [9,]  8 0.000007 0.000003
[10,]  9 0.000001 0.000000
[11,] 10 0.000000 0.000000
[12,] 11 0.000000 0.000000
[13,] 12 0.000000 0.000000
[14,] 13 0.000000 0.000000
[15,] 14 0.000000 0.000000
[16,] 15 0.000000 0.000000
[17,] 16 0.000000 0.000000
[18,] 17 0.000000 0.000000
[19,] 18 0.000000 0.000000
[20,] 19 0.000000 0.000000
[21,] 20 0.000000 0.000000

Above a table of this probabilities. To detect from observed data which of this two distributions one have, is what it takes to decide, in this case, if $N=100$ or if $N=20$. It is obviously quite hard, and the instability of the resulting estimators is only to be expected. This example also indicated that the instability is mainly for small $p$. You say you expect $p$ around 0.7, so the problem might be more stable then. You could investigate that for your data by finding the maximum likelihood estimator as a function of a known $p$, and plotting that for $p$ in some confidence interval. Or you could go full Bayes, this is a case where even some rather vague prior information could be helpful.

The parameters are indeed estimable. It is clear that $N_i \ge \max_t(x_{it})$, so it is possible to use that maximum count as an estimator of $N$. That estimator will be strongly consistent, and a parameter with a consistent estimator must be estimable. But, as the above example shows, the estimability is almost a formality; in practice distributions with very different $N$ are very close, so $N$ is very weakly estimable.

I will not give details of the estimation methods here, but give a few references you can check out:

Ingram Olkin, A John Petkau, James V Zidek: A comparison of N estimators for the Binomial Distribution. JASA 1981. This is a classic paper which develops and analyzes ML and moment estimators, and some stabler variants. It also shows, interestingly, that in many cases the method-of-moments estimator is better than the ML estimator!

Raymond J Carrol and F Lombard: A note on N estimators for the binomial distribution. JASA 1985.
Develops an alternative, stabler & maybe better estimator, based in integrating $p$ out of the likelihood. Also notes the lack of sufficiency of the summed counts.

J Andrew Royle: N_Mixture Models for Estimating Population Size from Spatially Replicated Counts. Biometrics, 2004. This gives another, alternative Bayesian approach which you may try.

Back to your concrete question. You SHOULD NOT sum the counts over your two regions! That will lose information. If you introduce $N=N_1 + N_2$ then the log-likelihood function can be written as a function of $N$, $p$ and $N_1$ (or $N_2$). Then the extra parameter $N_1$ should be eliminated by some procedure. I will come back to that, but no there is no time!

  • 1
    $\begingroup$ Thanks for your answer, unfortunately I cannot upvote. I'm sorry if I was unclear, but for the sum I mean the sum over different zones. I made an upgrade in my Q, adding an example (theoretical) and some maybe more understandable questions. $\endgroup$ Commented Nov 12, 2014 at 22:54
  • 1
    $\begingroup$ Despite this answer is fully detailed, I have still one doubt: suppose that I have one camera for each region,and suppose that the cameras are in neighbouring (but without flow) regions. I'm only interested in $N$ and not in the singles values $N_1$ and $N_2$. Is there a difference by taking only one big camera? The observed variables $y_{t1}=x_{t1}+x_{t2}$ of the big camera are still binomials with parameter $N,p.$ What I loose (it seem to me), is only information on the local behaviour, but this is not important for me. $\endgroup$ Commented Nov 13, 2014 at 22:02
  • 4
    $\begingroup$ That is not correct! Information about $N$ is contained in the variation in the individual counts! $\endgroup$ Commented Nov 13, 2014 at 22:41

I stumbled upon this somewhat old post by chance, and I was (at first) floored by @kjetilbhalvorsen's excellent answer.

It was surprising to me initially that by summing you should lose information about the sum $N_1+N_2$. This answer seeks to demonstrate that fact and provide a (preliminary) framework of numerical study on the difference between the two approaches.

The Result

Let's start with the good stuff. Following the original question, suppose that $x_1, \ldots, x_k \sim \text{Binom}(N_1, p)$ and $y_1, \ldots, y_k \sim \text{Binom}(N_2, p)$. The parameter of interest is $N= N_1+N_2$.

Naturally, since $z_i = x_i + y_i \sim \text{Binom}(N_1 + N_2, p)$, one approach is to apply a standard estimator to the $z$ sample to estimate $N$ as $$\hat N_z \stackrel{\scriptsize \triangle}{=}\hat N_\text{s:mme}(z_1, \ldots, z_k),$$ where $\hat N_\text{s:mme}$ is a stabilized method of moments estimator described below, but we could apply this for any estimator of $N$. But as the other answer points out, this loses information. Thus a secondary approach is to first estimate $\hat N_1$ and $\hat N_2$ using the same procedure, setting $$\hat N_{xy} \stackrel{\scriptsize \triangle}{=}\hat N_\text{s:mme}(x_1, \ldots, x_k) + N_\text{s:mme}(y_1, \ldots, y_k).$$

According to @kjetilbhalvorsen, the second estimator should perform better, and indeed it does -- at least for the simulation settings I chose ($N_1=200$, $N_2=500$, $p=0.05$, $k=30$). Note that both estimators seem to severely underestimate the true value (this could surely be improved by increasing $k$).

Comparison of Estimators

Review of Olkin et al. (1981)

Note that one could extend this analysis to maximum likelihood (MLE) or other estimators. But the Olkin paper shows that, especially for small samples (values of $k$), the method of moments (MME) approach usually outperforms MLE. So to keep things simple, we stick with a (stabilized) MME approach here. Briefly, the MME's of $N$ and $p$ are given by \begin{align*} \hat N_\text{mme} &= \frac{\hat\mu^2}{\hat\mu - \hat\sigma^2} \\ \hat p_\text{mme} &= \frac{\hat\mu}{\hat N_\text{mme}}, \end{align*} where $\hat\mu$ and $\hat\sigma^2$ are given (for now) by the sample mean and sample variance. This simple estimator has some obvious problems; when $\hat\mu < \hat\sigma^2$ we obtain a (clearly absurd) negative estimate for $N$. Worse still, if a sample contains all zero's (which is reasonable in many cases) then we have a "divide by zero" problem. To account for this, Olkin et al. propose the following procedure for stabilizing the MME estimators. A very small number $\epsilon$ is added to each of the successes in the sample. Recommendations for choosing $\epsilon$ are given in the paper, but I will not give details here. Similar adjustments are made to ensure that $\hat N$ is not less than any of the data points. In brief, the stabilized estimator for $N$ can be stated as follows \begin{align*} \phi &= \begin{cases} \hat\mu/\hat\sigma^2, & \hat\mu/\hat\sigma^2 \geq 1 + 2^{-1/2} \\[1.2ex] \text{max}\left\{1+\sqrt{2}, \left(x_{(k)} - \hat\mu\right)/\hat\sigma^2\right\}, & \hat\mu/\hat\sigma^2 < 1 + 2^{-1/2} \end{cases} \\ \hat N_\text{s:mme} &= \text{max}\left\{x_{(k)}, \frac{\hat\sigma^2\phi^2}{\phi-1}\right\} \end{align*} where $x_{(k)}$ is the largest of $x_1, \ldots, x_k$. The estimator (for all estimators $\hat N$) for $p$ is simply given by $\hat p = \hat\mu(z_1, \ldots, z_k)/\hat N$.

R Code

For anybody interested in extending the analysis, my code is given here. The function to produce $\hat N_\text{s:mme}$ is

N_mme_est <- function(xx){
  mu_hat <- mean(xx)
  s2_hat <- var(xx)
  if(mu_hat/s2_hat >= 1 + 1/sqrt(2)){
    phi <- mu_hat/s2_hat
    zz <- (max(xx) - mu_hat)/sqrt(s2_hat)
    phi <- max(z/sqrt(s2_hat), 1 + sqrt(2))
  N_hat <- max(max(xx), s2_hat*phi^2/(phi-1))

The following code runs a simple simulation study to compare the two cases.

N1 <- 200
N2 <- 500
N <- N1 + N2
p <- 0.05

M <- 5000 
N_hat_z <- N_hat_xy <- rep(NA, M)
for(i in 1:M){
  # Simulate data
  k <- 30 # Number of samples
  x <- rbinom(k, N1, p)
  y <- rbinom(k, N2, p)
  z <- x + y
  # Estimator based on z
  N_hat_z[i] <- N_mme_est(z)
  # Estimator based on x and y
  N_hat_x <- N_mme_est(x)
  N_hat_y <- N_mme_est(y)
  N_hat_xy[i] <- N_hat_x + N_hat_y

I used the following code to produce the figure above.

     col=adjustcolor("orange", alpha.f=0.4), 
     xlim=range(c(N_hat_z, N_hat_xy, N)),
     main="Comparison of Estimators",
     xlab="Value of N")
     col=adjustcolor("dodgerblue", alpha.f=0.4), 
     xlim=range(c(N_hat_z, N_hat_xy, N)),
abline(v=N, lwd=3, lty=1, col='firebrick')
legend("topright", c("Estimate based on z", "Estimate based on x & y", "True value"),
       fill=c(adjustcolor("orange", alpha.f=0.4),
              adjustcolor("dodgerblue", alpha.f=0.4),


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.