Formula for simulating radioactive decay for a large number of isotopes

Question

Currently I'm working on a project where I need to simulate the decay of a number of isotopes after each second.

One way to do so is each second do a uniform random roll for each particle, and if it is smaller than the decay constant, then count it as a decay.

Although this model is accurate it is very inefficient for large numbers of particles.

My second approach was to use a binomial distribution to model the decay. This works fine when the number of particles is sufficiently large, but fails when the probability of a decay is small.

For instance lets take Tritium, with a half life of $3.88524\times10^{8}$ seconds. If we have $2,000,000$ particles of Tritium, we expect 1 decay about every 5 minutes. By running the first simulation, indeed this is the case; however for the second one a decay never happens (for as long as I ran the experiment).

Any advice on how to better model this phenomenon?

EDIT: here is the first code that simulates correctly

production = 0;
for(var i = 0; i < number; i++){
    rand = Math.random();
    if(rand <  decay_constant){
        production++;
    }
}

Here is the code that simulates incorrectly

// p is the decay constant
var p = Math.log(2) / half_life;
var q = 1-p;
var mean = number*p;
var variance = number*p*q;
var std = Math.sqrt(variance);          
prod = Math.round(numberGenerator.nextGaussian()*std+mean);

Answer 1

Here are some exact probability computations comparing Poisson, binomial, and normal distributions. I believe the Poisson distribution is best as long as the period of time considered is several orders of magnitude less than the half life. (For longer periods of time adjustments of the kind suggested in the Comment by @A.S. are helpful.)

Poisson model for one minute with binomial approximation. The number of decays $X$ observed in one minute should average $\lambda = 1/5.$ So we use $X \sim Pois(\lambda = 1/5).$ Because the number of particles is $n = 2 \times 10^6$ a binomial distribution with $p = \lambda/n = 10^{-7}$ has the same mean. For these parameters the fit is quite good. The following R code makes a table of the respective probabilities for $x = 0, 1, 2, 3$ and the left panel of the figure at the end (Poisson bars, binomial open circles).

Roughly speaking, a sufficient condition for a reasonable fit between binomial and Poisson is that $n$ be very large, $p$ very small, and $\lambda = np$ of moderate size. However, this Poisson distribution is so severely skewed that it does not make sense to try a normal approximation.

 # ONE MINUTE
 x = 0:3;  lam = 1/5;  n = 2*10^6;  p = lam/n
 pdf.pois = dpois(x, lam);  pdf.bino = dbinom(x, n, p)
 cbind(x, pdf.pois, pdf.bino)
 ##     x    pdf.pois    pdf.bino
 ##     0 0.818730753 0.818730745
 ##     1 0.163746151 0.163746165
 ##     2 0.016374615 0.016374610
 ##     3 0.001091641 0.001091640

 par(mfrow = c(1,2))  # permits two panels per figure
 plot(x, pdf.pois, type="h", lwd=2, ylab="PDF", main="POIS(.2) With BINO Aprx")
 points(x, pdf.bino, col="blue")
 abline(h=0, col="green2")

Poisson model for one hour with binomial and normal approximations. The number of decays $X$ in one hour has a mean 60 times as large, so we use $X \sim Pois(\lambda = 12).$ The corresponding binomial has the same $n$ as above and $p = \lambda/n = 6 \times 10^{-6},\,$ so that the means match.

As $\lambda$ increases, Poisson distributions become less skewed and a normal approximation becomes feasible. You can see, by the plot in the right hand panel below, that $\lambda = 12$ is not yet quite large enough for a really good approximation. The appropriate normal approximation is $Norm(\mu = \lambda,\, \sigma =\sqrt{\lambda}).$ The additional R code below finishes the plot.

 # ONE HOUR   
 x=0:30;  lam=12;  m=2*10^6;  p = lam/n
 pdf.pois = dpois(x, lam);  pdf.bino = dbinom(x, n, p)
 plot(x, pdf.pois, type="h", lwd=2, ylab="PDF",main="POIS(12) W/ BINO & NORM Aprx")
 points(x, pdf.bino, col="blue")
 curve(dnorm(x, lam, sqrt(lam)), col="maroon", add=T) 
 abline(h=0, col="green2")
 par(mfrow = c(1,1))  # returns to default single-panel plots

As an example: for the probability of 10 or fewer events in an hour, the Poisson and binomial distributions both give 0.3472, and the normal approximation 0.3325 (all rounded to four places).

 ppois(10, 12);  pbinom(10, 2*10^6, 6*10^-6);  pnorm(10.5, 12, sqrt(12))  
 ## 0.3472294
 ## 0.3472288
 ## 0.3325028

---

Note: I hope the R code is sufficiently transparent that you can follow it as much as you need to, and perhaps make minor modifications. R software is available free of charge and with easy installation at www.R-project.org. There are Windows, Mac, and UNIX versions. I have used only the 'base' of R (no libraries) and put #-symbols to keep answers and comments from interfering with running the code as shown.