Why is my exponential decay not working?

Question

I'm simulating exponential decay. Why is it not working? Here's the code

from math import *
from random import random
from time import sleep

class atom():
    def __init__(self, pos):
        self.pos = pos
        self.isalive = True

class substance():
    def __init__(self, halflife):
        self.halflife = float(halflife)
        self.rate = float(log(2) / self.halflife)
        self.life = float(0)
        self.probdecay = float(0)
        self.struct = [atom(i) for i in range(100)]
        self.isalive = True
        print self.rate

    def sim(self):
        print self.probdecay
        for i in range(len(self.struct)):
            if self.struct[i].isalive:
                if random() < self.probdecay:
                    self.struct[i].isalive = False

        self.life += 1.0
        self.probdecay = 1.0 - pow(1.0 - self.rate, self.life)

    def getstr(self):
        s = ""
        for i in self.struct:
            if i.isalive:
                s += "o"
            else:
                s += " "
        return s

s = substance(5)

while s.isalive:
    s.sim()
    print s.getstr()
    sleep(1)

Here's the output

0.138629436112
0.0
oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
0.138629436112
oooooo o ooooooooooo oooooooo ooooooo o oooooooooooo oooooo oooooooo oooo ooo oooooooo  oooo   oooo 
0.258040751667
o  ooo o   oo oo ooo o    ooo o ooooo o oo  oooooo o   oooo  oo oooo oo   ooo ooo  ooo   ooo   oo   
0.360898143882
   o o     oo     oo o    ooo   ooo o   oo  o o oo     o o    o    o o    ooo o o   oo   o o   oo   
0.449496473814
   o        o        o    o o    o      o     o        o      o      o     o                   oo   
0.525812467226
   o        o                    o      o     o                      o                              
0.591548817506
            o                           o                                                           
0.648172174615

0.696945867656

0.738958091134

0.775146183762

Obviously, that is not exponential decay with halflife of 5. What's going on here?

Answer 1

It looks like you mixed up the probability of a particle decaying on a certain step, and the probability of a particle decaying by a given step. The probability of a particle decaying on any given step is the same, regardless of which step it is, while the probability that a particle has decayed by a given step is monotonically increasing. Instead of gradually increasing probdecay, set it once, when you initialize the substance:

def __init__(self, halflife):
    self.probdecay = 1 - 0.5 ** (1.0 / halflife)
    ...