8
$\begingroup$

Situation:

Say I have a Poisson process, like radioactive decay, producing R particles per second. I measure with a detector. There is a probability P that a particle will be detected by the detector.

Things I think I know:

  1. The inter-arrival time of the particle emission is exponentially distributed with parameters based on R.
  2. The number of particles emitted before detection is given by a negative binomial based on P.
  3. If a number N is sampled from (2), a single sample of inter-arrival time for detected particles can be given by the sum of N samples from (1). This sum can be obtained by sampling from a gamma distribution with parameters based on N and R.

My question:

If a single inter-arrival time can be calculated by sampling from a gamma based on N and R, how does the number of detector counts in an interval end up being Poisson again? (To be Poisson, the inter-arrival time for the detector must be exponential, not distributed according to some weird gamma thing.) Of course N is fluctuating, but I can't see how this works out.

However, I am almost completely sure detector counts are in fact Poisson distributed. Could somebody show me the math? Thanks for the help!

EDIT:

I found this paper: Fried, D. L. "Noise in photoemission current." Applied Optics 4.1 (1965): 79-80.

Which shows the result that a binomially selected poisson random variable is also Poisson with a rate given by PR. This confirms the comment by jbowman. Still, I would be interested in seeing the explanation of how my process of generating the inter-arrival time at the detector using the negative binomial and gamma distribution is incorrect. This is my major mental hiccup. Thank you.

EDIT 2:

I wrote this matlab script to test whether what I was trying with the gamma distribution worked. Turns out that somehow the gamma inter-arrival times generated with a geometrically distributed N are exponential and agree with the inter-arrival times suggested by Poisson(PR). (ia2 and ia3 are identically distributed). Any idea how this works out analytically? It was not intuitively obvious to me!

close all
n = 100000;
ia1 = exprnd(1,n,1); % create exponentially distributed inter-arrival times
t1 = cumsum(ia1); % running sum (the real experiment time)

mask = (rand(n,1) > 0.5); % flip a coin
t2 = t1(mask); % get only the events for which "the coin landed on heads"
ia2 = diff(t2); % calculate the inter-arrival times at the detector.

% plot the distributions
figure; hist(ia1,100); title('exponential inter-arrival times');
figure; hist(ia2,100); title('binomial sampled inter-arrival times');

%%
spacing = geornd(0.5,n,1) + 1; % how many events before we get heads
ia3 = gamrnd(spacing,ones(n,1)); % generate the interarrival times with gamma
figure; hist(ia3,100); title('geom/gamma inter-arrival times');
$\endgroup$
  • $\begingroup$ #2 is actually not correct; if each particle has a probability $P$ of being detected, the distribution of detected particles per second is Poisson($RP$) (assuming detections are independent etc.) For #3, if $N$ is sampled from 2, then you don't have a single sample of inter-arrival times; you have a single observation of the sum of $N$ interarrival times, which is indeed distributed Gamma with shape parameter $N$. Consequently, the premise of your question ("If a single interarrival time...") isn't true. $\endgroup$ – jbowman May 29 '14 at 17:28
  • $\begingroup$ I don't understand how you know the rate is Poisson(RP). Could you show me? That is the very heart of this question I think. In #2, I suppose that I have if I have a P chance to hit the detector, the number of particles emitted before hitting the detector is Geometrically distributed with a mean of 1/P. Thus, I can calculate can sample from this Geometric distribution to get N, then sum up N inter-arrival times to get a single inter-arrival time at the detector. Can you explain the flaw in this logic? I think your statement about the rate being Poisson(RP) is important. Thank you! $\endgroup$ – user487100 May 29 '14 at 18:14
  • $\begingroup$ Are you somewhat familiar with moment generating / characteristic functions? I'd write it out using that approach, as it's simple, unless it's also unhelpful. $\endgroup$ – jbowman May 29 '14 at 18:21
  • $\begingroup$ No I have not worked with moment generating functions. Do you have some idea of how to show that Poisson + some fixed probability of acceptance just scales the poisson rate? I am willing to learn the moment generating function approach based if you could show how this works out. $\endgroup$ – user487100 May 29 '14 at 18:43
  • 1
    $\begingroup$ It'll be much later today (Pacific Std Time), I'm afraid; I can do it the straightforward way too, which will be less opaque. $\endgroup$ – jbowman May 29 '14 at 19:57
4
$\begingroup$

A quick non-technical argument might use Jackson networks. In your case total external arrivals is rate $R$, and there are no internal transitions (observed particles don't switch to the unobserved queue). The splitting proportion between the observed and unobserved nodes $p_{0i}$ is $P$, so the

$\lambda_{obs}=RP$

If you're looking for first principles, call $O(t)$ the observed counting process, and $N(t)\sim PP(r)$ the total counting process. Where each arrival in $N(t)$ gets logged in $O(t)$ with probability $p$. So that if for some $s$ we have $N(s)=n$ then $O(s)$ has a binomial($n,p$) distribution.

This approach uses probability generating functions:

$E[z^{O(t)}|N(t)=n]=\sum_{j=0}^{n}z^{j} {n \choose j}p^j(1-p)^{n-j}=(1-p+pz)^{n}$

Last equality by the binomial theorem. Then, unconditionally, since $N(t)\sim Poisson(rt)$:

$E[z^{O(t)}]=E[E[z^{O(t)}|N(t)=n]]=\sum_{n=0}^{\infty}(1-p+pz)^{n}\frac{rt^n}{n!}e^{-rt}=e^{-rt}e^{rt(1-p+pz)}=e^{rpt(z-1)}$

Which is the probability generating function of a Poisson($rpt$) random variable.

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