1
$\begingroup$

I'm a lowly physicist, so I hope you will forgive me if I botch some terminology and notation here.

I perform an experiment where on each trial I start with some number of particles $N$ that is approximately Poisson distributed, say $N \sim {\rm Pois}(\overline{N})$. I then perform some action that bins these particles into two groups, $a$ and $b$, such that $N_a \sim B(N,p)$ and $N_b = N-N_a$. The value of $p$ is unknown to me. I know how to calculate ${\rm Var}[N_a]$, ${\rm Var}[N_b]$, and ${\rm Cov}[N_a,N_b]$ using the laws of total variance and covariance. I find $${\rm Var}[N_a] = \overline{N} \, p(1-p) + p^2 {\rm Var}[N],$$ $${\rm Var}[N_b] = \overline{N} \, p(1-p) + (1-p)^2 {\rm Var}[N],$$ $$ {\rm Cov}[N_a,N_b] = -\overline{N} \, p(1-p) + p(1-p) {\rm Var}[N]. $$

Based on measurements of $N_a$ and $N_b$, I would like to estimate the quantity $Y=2p-1$. I would like to consider two scenarios; one in which I can measure both $N_a$ and $N_b$ in a single experiment (i.e., for a single value of $N$), and another in which I cannot. My question pertains to the choice of estimator for $Y$ in the second scenario.

  1. If I can measure both $N_a$ and $N_b$ in a single trial, I can estimate $Y$ by $$ \hat{Y} = \frac{N_a-N_b}{N_a+N_b}, $$ and approximating the variance in $\hat{Y}$ by a series expansion, I find $$ {\rm Var}[\hat{Y}] = 4 \frac{N_b^2 {\rm Var}[N_a] + N_a^2 {\rm Var}[N_b] - 2 N_a N_b {\rm Cov}[N_a,N_b]}{(N_a+N_b)^4} $$ which gives $$ {\rm Var}[\hat{Y}] = \frac{4p(1-p)}{\overline{N}}. $$ As I understand from this SE post, this is the expected variance for the MLE of $\hat{Y}$. This estimator of $\hat{Y}$ normalizes out fluctuations in $N$; nowhere does ${\rm Var}[N]$ appear in this expression. This is desirable, from my point of view.

  2. What if I can only measure one of $N_a$ or $N_b$ in a single trial of my experiment? I can imagine several ways of estimating $Y$, for example I could simply use the same estimator as above, but if $N_a$ and $N_b$ do not derive from the same sample of $N$, they are uncorrelated and ${\rm Var}[N]$ no longer drops out of ${\rm Var}[\hat{Y}]$. One reasonable possibility seems like $$ \hat{Y}_0 = \frac{N_a - N_b}{\overline{N}}, $$ where I estimate $\overline{N} = \overline{N_a}+\overline{N_b}$ from a sample mean over several trials. I can also dream up something goofy like $$ \hat{Y}_1 = \frac{N_a}{N_a + \overline{N_b}} - \frac{N_b}{\overline{N_a}+N_b}, $$ which from some MC simulations I've performed, seems to be biased, but it's not clear to me why.

My question is, in scenario 2, which estimator for $Y$ should I choose, and why? I would like to estimate $Y$ with small bias and with the minimum variance possible.

Edit: I do not know $\overline{N}$ ($=\lambda$); I can only measure $N_a$ and $N_b$. If my best estimator for $p$ is $N_a/\lambda$, how should I estimate $\lambda$?

$\endgroup$
  • $\begingroup$ How many trials will you perform? just one? Under scenario 2, could you get the sum of N if multiple trials will be performed? $\endgroup$ – user158565 Jun 8 '19 at 2:30
  • $\begingroup$ I can perform many trials, and in scenario 2 I typically measure N_a on one trial, and N_b on the next. I do some large number of trials, and I suppose my best estimate for lambda is N_a + N_b, averaged over trials. $\endgroup$ – Will C Jun 8 '19 at 3:23
  • $\begingroup$ Google 'splitting' and 'thinning' Poisson processes as here. $\endgroup$ – BruceET Jun 8 '19 at 11:52
1
$\begingroup$

Suppose $N \sim \mathsf{Pois}(\lambda).$ (I'm using $\bar N = \lambda.)$ According to your description it seems you observe $N_a \sim \mathsf{Pois}(p\lambda).$ Then $E(N_a) = p\lambda,$ so $\hat p = N_a/\lambda$ is the natural estimator $p.$

One reasonable 95% confidence interval for $p\lambda$ is $$N_a + 2 \pm 1.96\sqrt{N_a + 1}.$$ Dividing by $\lambda$ we get the endpoints of the above 95% CI for $p.$

The CI for $p\lambda$ is based on a normal approximation, which works reasonably well if $N_a \ge 50.$

For example if $\lambda = 500$ and $N_a = 288$ then the CI for $p\lambda$ is $$290 \pm 1.96\sqrt{289}.$$ Dividing by $500$ we get $(0.513, 0.647)$ as a 95% CI for $p$ and point estimate $\hat p \approx 0.58.$

If events are rare so that you have much smaller $N_a,$ then you need to look into a CI for $p\lambda$ that does not use a normal approximation.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ I'm not sure I understand the question, but I'll try to answer: I create an unknown number of particles N on each trial. I then "turn them" into a-type with probability p. All that are left are b-type. The problem is that I can't measure both N_a and N_b on a given trial. $\endgroup$ – Will C Jun 8 '19 at 3:25
  • $\begingroup$ The typical N, although I can't measure it directly, is around 1000. The probability p actual varies as a function of other parameters anywhere between 0 and 1. $\endgroup$ – Will C Jun 8 '19 at 3:41
  • $\begingroup$ Thanks for additional info. I simulated several models, all with results consistent with my Answer. I'm taking $\lambda = \bar N$ as a constant, $N$ and $N_a$ as random variables. Will look at it again tomorrow. $\endgroup$ – BruceET Jun 8 '19 at 4:25
1
$\begingroup$

Scenario 2:

Suppose you perform $n = n_1 + n_2$ trials and observe $N_a$ in $n_1$ trials and $N_b$ in $n_2$ trials. No trial has both $N_a$ and $N_b$. Trails are independent.

Let $N_A = \sum N_a$ and $N_B=\sum N_b$. Then $N_A \sim Poisson(n_1p\bar N)$ and $N_B \sim Poisson(n_2(1-p)\bar N)$

Then you can right the joint log likelihood of $N_A$ and $N_B$, and get the $\hat p$ MLE of $p$ and its CI. It is straightforward to get $\hat Y$ and its CI, because $Y$ is linear function of $p$.

| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.