## Stochastic model simulation: radioactive decay as
## the 1st event in a Homogeneous Poisson Process for...
N0 = 5000 # particles
## For each, sample an exponentially distributed
## random number with rate r
r = 1
## Save those event times in a vector (and sort them)
Ts = sort(rexp(N0,r))
## Show number not decayed as a function of time...
N = c() # number at time t not-yet-decayed.
Time = seq(0, max(Ts), length=500)
N[1]=N0
for(i in 2:length(Time)) {
N[i] = sum(Ts > Time[i]);
}
plot(Time, N)
## Add the "mean field" model exponential decay curve
curve(N0*exp(-r*x),0,max(Time),col="orange", lwd=2, add=TRUE)