Decay chain simulated in Python

Question

I am trying to make a computer program which can model the decay of 1,000,000 atoms of Radon-222 (into Polonium-218 into Lead-214 into Bismuth-214 into Lead-210). Each atom has a 0.000126 chance of decaying into Polonium-218 after one unit time. I have a list of 1,000,000 0's m assigned to the variable Rn to represent the Radon-222. Then, to decide if a particular element decays, I choose a random integer between 1 and 1,000,000: if this number is less than 126, the value of Rn[j] switches from 0 to 1 (decay), and we add a 0 to the list Po for polonium, etc. So after 1 unit time, I should expect there to be around 126 Polonium atoms, after another unit time, we should have around 126 extra Polonium atoms, but some of these could have decayed into Lead-214.

We treat time as a discrete entity in this modelling. See here for a list of ideal values after 10 unit times.

Here is a code which works, but is rather slow.

import random

def one(s):
    count=0
    for i in range(len(s)):
        if s[i]==1:
            count+=1
    return count

Rn=[1]*1000000
Po=[0]*1000000
Pb214=[0]*1000000
Bi=[0]*1000000
Pb210=[0]*1000000

print([0,1000000,0,0,0,0])

for i in range(1,100):
    for j in range(len(Rn)):
        rndecay=random.random()
        if rndecay<=0.000126:
            Rn[j]=2        

    for j in range(len(Po)):
        if Po[j]==1:
            podecay=random.random()
            if podecay<=0.203:
                Po[j]=2

        if Rn[j]==2 and Po[j]==0:
            Po[j]=1

    for j in range(len(Pb214)):            
        if Pb214[j]==1:
            decay214=random.random()
            if decay214<=0.0255:
                Pb214[j]=2
        if Po[j]==2 and Pb214[j]==0:
            Pb214[j]=1

    for j in range(len(Bi)):
        if Bi[j]==1:
            bidecay=random.random()
            if bidecay<=0.0346:
                Bi[j]=2
        if Pb214[j]==2 and Bi[j]==0:
            Bi[j]=1

    for j in range(len(Pb210)):
        if Bi[j]==2 and Pb210[j]==0:
            Pb210[j]=1

    print([i,one(Rn),one(Po),one(Pb214),one(Bi),one(Pb210)])

How can I optimize this code?

Answer 1

I think you might have better luck approaching the problem differently. Instead of having 1,000,000 element lists, could you just keep track of the number of atoms of each species?

The number of decays in a population of n with a decay probability p follows a binomial distribution with those parameters. Many basic statistics packages will allow you to generate random numbers that follow a distribution (rather than the uniform distribution used by random.random). Thus, you only need to pull one random value for each possible transition and then update the counts with the results. This methodology is exactly equivalent to the one employed in the original post.

Here's an example using counts with transition probabilities that I made up:

from numpy.random import binomial

atoms = {'Rn222': 1000000,
         'Po218': 0,
         'Pb214': 0,
         'Bi214': 0,
         'Pb210': 0}

for _ in range(100):
    Rn222_Po218 = binomial(atoms['Rn222'], 0.000126)
    Po218_Pb214 = binomial(atoms['Po218'], 0.001240)
    Pb214_Bi214 = binomial(atoms['Pb214'], 0.003450)
    Bi214_Pb210 = binomial(atoms['Bi214'], 0.012046)

    atoms['Rn222'] -= Rn222_Po218
    atoms['Po218'] += Rn222_Po218

    atoms['Po218'] -= Po218_Pb214
    atoms['Pb214'] += Po218_Pb214

    atoms['Pb214'] -= Pb214_Bi214
    atoms['Bi214'] += Pb214_Bi214

    atoms['Bi214'] -= Bi214_Pb210
    atoms['Pb210'] += Bi214_Pb210

print atoms

Edited to add example. Edited to add binomial explanation from comment.