Poisson Process – Inter-arrival Time Distributions¶

Principle¶

A homogeneous Poisson process with rate λ can be simulated by:

  1. Generate N ≈ λ·T independent uniform random numbers on [0,T]
  2. Sort them → arrival times t₀ < t₁ < … < t_{N-1}

The time difference between events separated by k steps
Δtₖ = tᵢ – t_{i-k}
follows a Gamma(k, λ) distribution (Erlang when k integer):

$p(Δtₖ) = \frac{λ^k · (Δt)^{k-1} · e^{-λ Δt}}{(k-1)!}$

In [2]:
// 1 — Generate Poisson events 
#include <vector>
#include <algorithm>
#include <TRandom3.h>
#include <cstdio>
using namespace std;

vector<double> gTimes;           // global vector — automatically shared between cells
gTimes.clear();                  // make sure it's empty

double rate    = 100.0;          // events per second
double total_T = 10000.0;        // observation time
double xmax    = 0.20;           // x-axis limit: 0 to 0.2 s

gRandom = new TRandom3(0);
Long64_t N = (Long64_t)(rate * total_T + 0.5);

gTimes.reserve(N);
for(Long64_t i = 0; i < N; i++)
    gTimes.push_back(gRandom->Uniform(total_T));

sort(gTimes.begin(), gTimes.end());

printf("Generated %lld events → gTimes.size() = %zu\n", N, gTimes.size());

Generated 1000000 events → gTimes.size() = 1000000
In [3]:
//  Plot inter-arrival distributions
double xmax = 0.20;
int    bins = 800;
double bw   = xmax / bins;
Long64_t N  = gTimes.size();

TH1D* h[5];
int color[5] = {kRed+1, kGreen+2, kBlue+1, kBlack, kMagenta+2};

for(int k = 0; k < 5; k++)
    h[k] = new TH1D(Form("h%d",k+1), "", bins, 0, xmax);

for(size_t i = 0; i < gTimes.size(); i++)
    for(int k = 1; k <= 5; k++)
        if(i >= (size_t)k)
            h[k-1]->Fill(gTimes[i] - gTimes[i-k], 1.0/N/bw);

TCanvas c("c", "Poisson Inter-arrival Times (0–0.2 s)", 860, 560);
c.SetLogy(); c.SetLeftMargin(0.12);

for(int k = 0; k < 5; k++) {
    h[k]->SetLineColor(color[k]);
    h[k]->SetLineWidth(2);
    h[k]->Draw(k==0 ? "" : "SAME");
}

TF1* f = new TF1("f", "[0]*exp(-[1]*x)", 0, xmax);
f->SetParameters(100, 100);
f->SetLineColor(kOrange+7);
f->SetLineStyle(2);
h[0]->Fit(f, "R");

h[0]->SetTitle("");
h[0]->GetXaxis()->SetTitle("Inter-arrival time Δt (s)");
h[0]->GetYaxis()->SetTitle("Probability density (s^{-1})");

gStyle->SetOptStat(0);
gStyle->SetOptFit(1111);
c.Draw();

****************************************
Minimizer is Minuit2 / Migrad
Chi2                      =       409.66
NDf                       =          426
Edm                       =  2.25653e-06
NCalls                    =           35
p0                        =      100.097   +/-   0.141315    
p1                        =      100.132   +/-   0.0997951