# p is the probability that a particle will decay on a give time step
p=0.05
# sample_size is the number of particles we will simulate the decay of
# a value of 1 will correspond to undecayed and 0 will correspond to decayed
sample_size=100
# decay_time is how many steps in the simulation
decay_time=50
# start out with all particles undecayed 
# -- remeber there are sample-size particles
particles=rep(1,sample_size)
# making a place to hold the sum of undecayed particles after each step
# -- remember each particle will take decay_time steps
number_left=integer(decay_time)

# go through each step of time 
# -- known as a loop
for(i in 1:decay_time){
  # use runif (randon-uniform) to obtain an array of pseudo-random between 0 and 1 
  rand=runif(sample_size,0,1)
  # if the random number for that particles is less than p multiple it 
  # by 0 and turn it into a decayed particles
  # that is, multiply particles by ifelse(rand<p, 0, 1) and make it 
  # the new particles
  particles = particles*ifelse(rand<p, 0, 1)
  # record the sum of the particles -- since the decayed particles are 
  # 0 the sum gives the number of undecayed particles 
  number_left[i]=sum(particles)
}
# plot the number of particles  -- supposed to "decay" exponentially
plot(number_left)