-- Generate Poisson random inter-event times.
-- Since the Poisson distribution is a simple exponential, the inverse of its integral can be used to transform the uniform random distribution.

poissonrand(aveper) = -ln(rand(1))*aveper

-- Check the distribution by accumulating a histogram.
include "::incl:histogram" -- Or use "histogram" XFun

-- make an array of random inter-event times
nevents = 600
aveper = 2.5
dt[i] = poissonrand(aveper) dim[nevents]

-- accumulate a histogram of number of events vs. inter-event time.
nbins = 40
histo(dt,Xmin,Xmax,nbins):
--histogram(dt,Xmin,Xmax,nbins) --XFun

plot bins
label mean:2.377
label aveper:2.500
Xmin=0; Xmax=5*aveper
Ymin=0;

-- plot the continuous Poisson distribution
poisson(t) = exp(-t/aveper)

plot bin1*poisson(X)

bint = (Xmax-Xmin)/nbins
bin1 = nevents*(1-poisson(bint))
