Estimate sample size if two properties are measured and their probabilities are known

Question

I'm sorry if this is some kind of standard exercise, but I just don't know how to search for it.

I have infinite amount of samples and I can measure two properties $A$ and $B$ for each sample. Chance to measure $A$ is $P(A)$ and for $B$ it's $P(B)$, both properties are independent. I take N samples and on $a$ of those I measure property $A$, but property $B$ is encountered $b$ times.

If $P(A)$, $P(B)$, $a$ and $b$ are known, how can I get $N$ and what's the error/deviation for this quantity?

I can also formulate the question as a problem with some numbers: I have infinite amount of boxes, in each box there is a candy wrapped in paper. 34% of the candies are made of chocolate and 10% of the papers are blue.

I've opened some boxes and seen 16 blue papers and have eaten 45 chocolates. How many boxes have I opened?

What changes if $P(A)$ and $P(B)$ are not precise but experimentally measured from $N_A$ and $N_B$ measurements instead? That means - before my measurements my brother opened some $N_A=522399$ boxes and that's where the $P(A)=34%$ came from. He's colourblind tho, but my sister saw him opening some $N_B=92587$ boxes and noticed that 10% of the papers were blue. Any changes if there are more than two properties?

Answer 1

We will approach this problem using maximum likelihood estimation to estimate $N$.

First, let us write out the probability of seeing $x$ chocolates and $y$ blue papers, given $N$ boxes were opened, the probability of seeing a chocolate is $p_x$, and the probability of seeing a blue paper is $p_y$. This is just the product of two binomial distributions with the same $N$:

$P(x,y) = \frac{N!}{x!(N-x)!} p_x^x(1-p_x)^{N-x}\frac{N!}{y!(N-y)!}p_y^y(1-p_y)^{N-y}$

We solve for the $N$ which gives us the highest probability of seeing the observations we actually saw, in this case 45 chocolates and 16 blue papers. This is better done by maximising the log of $P(x,y)$ for reasons of stability and avoiding really large or small numbers.

Here's some R code to do this in a rather brute-force way:

log.pxy <- function(N, px, py, x, y)
{
  dbinom(x, N, px, log=TRUE) + dbinom(y, N, py, log=TRUE)
}

results <- rep(-Inf,200)
for (N in 50:200) results[N] <- log.pxy(N, 0.34, 0.1, 45, 16)
which.max(results)
[1] 137

Our maximum likelihood estimate of $N$ is 137.

As for getting a confidence interval, we can use the fact that -2 * the observed log of the likelihood function ($P(x,y)$ viewed as a function of $N$ instead of $(x,y)$) is asymptotically distributed $\chi^2(1)$. We find the $N$ for which the log likelihood is "too far" below the maximum, using the $\chi^2(1)$ as our guide, and constructing a 95% confidence interval:

> min(which(2*results > max(2*results)-qchisq(0.95,1)))
[1] 111
> max(which(2*results > max(2*results)-qchisq(0.95,1)))
[1] 168

So a 95% confidence interval would be [110, 169] - we need to expand the interval by 1 on either end to get "outside" the qchisq 95% range.

As for your other questions: If there are more than two properties, you can expand the solution methodology in the obvious way and it will still work.

The more complex situation is when $p_x$ and $p_y$ are themselves estimates based on samples. I'd jointly estimate $p_x$, $p_y$, and $N$ in that case, which can be done with a nested procedure that iterates over $N$ in an outer loop, and, in an inner optimization, estimates $p_x$ and $p_y$ given all the data and $N$. (As you will see, this is quite straightforward; given $N$, it's just the usual estimate of $p_x$ and $p_y$.) We then find the $N$ which maximizes the log likelihood as before.

The likelihood function is more complex. Let's denote the other information with $N_a$ and $N_b$ sample sizes and observed values $x_a$ and $y_b$. We have:

$L(N,p_x,p_y) = {{N}\choose{x}}p_x^x(1-p_x)^{N-x} {{N_a}\choose{x_a}}p_x^{x_a}(1-p_x)^{N_a-x_a} \dots$

where the $\dots$ save us from writing out the $p_y$ part. Obviously we can combine some terms, but this form makes it a little easier to see what's going on.

Now for the R code. I'll assume, for concreteness, that $N_a=200$ and $x_a=68$, giving the point estimate for $p_x=0.34$, and $N_b=100$, $y_b=10$, giving the point estimate for $p_y=0.1$.

log.ll <- function(px, py, N, x, y, Na, Nb, xa, yb) {
  dbinom(x, N, px, log=TRUE) + dbinom(y, N, py, log=TRUE) +
      dbinom(xa, Na, px, log=TRUE) + dbinom(yb, Nb, py, log=TRUE)
}

x = 45
y = 16
Na = 200
Nb = 100
xa = 68
yb = 10

log.ll.N <- rep(-Inf,200)
for (N in 51:200) {
  px.hat <- (x+xa)/(N+Na)
  py.hat <- (y+yb)/(N+Nb)

  log.ll.N[N] <- log.ll(px.hat, py.hat, N, x, y, Na, Nb, xa, yb)
}

And, for the answers:

> which.max(log.ll.N)
[1] 135
> min(which(2*log.ll.N > max(2*log.ll.N)-qchisq(0.95,1)))
[1] 103
> max(which(2*log.ll.N > max(2*log.ll.N)-qchisq(0.95,1)))
[1] 180
>

For a slightly different point estimate of 135 for $N$, and a wider confidence interval of 102 - 181, as befits our new lack of precision about $p_x$ and $p_y$. We can recover our new estimates of $p_x$ and $p_y$ based on our combined sample:

> N <- 135
> (x+xa)/(N+Na)
[1] 0.3373134
> (y+yb)/(N+Nb)
[1] 0.1106383

I should also point out that our confidence interval is based on the profile log likelihood, not the log likelihood, but it's still a perfectly valid confidence interval.

Answer 2

This reminds me of the mark and recapture strategy to estimate a population size. This is modelled by the hypergeometric distribution, but in your case each property is a binomial.

For one property you would write the likelihood of observing $a$ as

$${{N}\choose{a}} P(A)^a(1-P(A))^{N-a}.$$

Because the properties are independent you can write the whole likelihood as the product

$$L(N) = {{N}\choose{a}} P(A)^a(1-P(A))^{N-a} {{N}\choose{b}} P(B)^k(1-P(B))^{N-b},$$

which you have to maximize for $N$. I would solve this using R because taking the log and differentiating with respect to $N$ introduces the digamma function (the derivative of the log-Gamma function), which does not make matters simpler.

Here is an R script that would give the estimate.

p <- .34  # P(A)
q <- .1   # P(B)
a <- 45
b <- 16
# Minimum possible value for N.
N <- max(a,b)
old_loglik <- -Inf
while (TRUE) {
   loglik <- dbinom(a, size=N, prob=p, log=TRUE) +
       dbinom(b, size=N, prob=q, log=TRUE)
   if (old_loglik > loglik) {
      # Passed the optimum.
      N <- N-1
      break
   }
   else {
      old_loglik <- loglik
   }
   N <- N+1
}
print (N)
# Answer is 137.

If $P(A)$ and $P(B)$ are not known and you have to estimate their value, then the initial problem becomes meaningless because you need to know $N$ in order to estimate them (anybody correct me if you know an alternative estimator).

If they were estimated on a first sample of size $n$ and you are drawing a new sample of size $N$ then the problem is more involved because you have to write the likelihood of everything and maximize it, this time with 3 parameters $(p, q, N)$.