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

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?

• Did you mean independent or mutually exclusive. Mutually eclusive means that they cannot both happen on the same draw. Independent means that if A occurs it has no influence on whether or not B occurs. Independent events could both happen on the same draw. So if they are mutually exclusive in your numerical example you would know that N is at least 61. If they are independent and P(A) and P(B) are both >0 then they are not mutually exclusive and you only know that N is at least 45. Now N is a random variable. So what do you mean by "how can I get N"? – Michael R. Chernick May 28 '12 at 14:12
• Are you wanting to derive its distribution or just compute its mean value and standard deviation? – Michael R. Chernick May 28 '12 at 14:13
• I mean independent as colour of paper and the filling is independent from one another. I need only mean value and deviation, but it would be interesting to see if it's possible to derive distribution. – Juris May 28 '12 at 20:49
• So independent in the sense that P(A∩B)=P(A) P(B). – Michael R. Chernick May 28 '12 at 20:55

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)
}


> 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.

• Thanks! I edited the question to explain the aditional question better. Basically it's about case where the given $P(A)$ and $P(B)$ have their own uncertainty. – Juris May 28 '12 at 21:03
• Ah, OK. I should have seen that. – jbowman May 28 '12 at 22:24

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)

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)$.