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. I'm not sure how your other question works; if $P(A)$ is estimated from $N_A$, then you know $N_A$ so there's nothing left to estimate - I think I must not be understanding the question.

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.

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(results > max(results)-qchisq(0.95,1)))
[1] 101
> max(which(results > max(results)-qchisq(0.95,1)))
[1] 183

So a 95% confidence interval would be [100, 184] - 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. I'm not sure how your other question works; if $P(A)$ is estimated from $N_A$, then you know $N_A$ so there's nothing left to estimate - I think I must not be understanding the question.

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. I'm not sure how your other question works; if $P(A)$ is estimated from $N_A$, then you know $N_A$ so there's nothing left to estimate - I think I must not be understanding the question.