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:
- The inter-arrival time of the particle emission is exponentially distributed with parameters based on R.
- The number of particles emitted before detection is given by a negative binomial based on P.
- 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.
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!
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.
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');