|
| 1 | +# ## The basics of Point Process Modeling |
| 2 | +# |
| 3 | +# This tutorial provides a brief introduction to point process modeling using PointProcesses.jl. |
| 4 | +# We will cover the basic math of point processes, how to define basic models, simulate data, |
| 5 | +# and fit models to observed data. For further, and more advanced pedagogy, please refer to |
| 6 | +# any of the following great resources: |
| 7 | +# |
| 8 | +# - [DaleyJones2013](@cite) |
| 9 | +# - [Cox2018](@cite) |
| 10 | +# |
| 11 | +# The focus of this tutorial is intentionally narrow: we aim to build intuition by starting |
| 12 | +# with the simplest point process model and making the connection between |
| 13 | +# mathematical definitions, simulation, and likelihood-based inference explicit. |
| 14 | + |
| 15 | +# ### Basic Point Process Concepts |
| 16 | +# |
| 17 | +# A point process is a mathematical model for random events occurring in time or space. |
| 18 | +# In this tutorial, we focus on **temporal point processes**, where events are represented |
| 19 | +# as a set of event times: |
| 20 | +# |
| 21 | +# {t₁, t₂, …, tₙ} ⊂ [0, T] |
| 22 | +# |
| 23 | +# The key components of a point process are: |
| 24 | +# |
| 25 | +# - **Events**: The occurrences we are modeling, represented as points in time. |
| 26 | +# - **Counting Process**: N(t), the number of events that have occurred up to time t. |
| 27 | +# - **Intensity Function**: A (possibly time- and history-dependent) function λ(t | 𝓗ₜ) |
| 28 | +# describing the instantaneous rate at which events occur. |
| 29 | +# - **History**: 𝓗ₜ, the record of events strictly before time t. |
| 30 | + |
| 31 | +# Intuitively, the intensity function λ(t | 𝓗ₜ) satisfies |
| 32 | +# |
| 33 | +# λ(t | 𝓗ₜ) Δt ≈ P(event in [t, t + Δt) | 𝓗ₜ) |
| 34 | +# |
| 35 | +# for small Δt. Different choices of λ give rise to different classes of point processes: |
| 36 | +# |
| 37 | +# - **Homogeneous Poisson process**: λ(t) = λ (constant) |
| 38 | +# - **Inhomogeneous Poisson process**: λ(t) varies with time. |
| 39 | +# - **History-dependent processes**: λ(t | 𝓗ₜ) depends on past events (e.g. Hawkes processes) |
| 40 | + |
| 41 | +# ### The Homogeneous Poisson Process |
| 42 | +# |
| 43 | +# The simplest and most widely used point process is the **homogeneous Poisson process**. |
| 44 | +# It is defined by a single parameter λ > 0, the constant event rate, and is characterized by: |
| 45 | +# |
| 46 | +# 1. **Independent increments**: The numbers of events in disjoint time intervals are independent. |
| 47 | +# 2. **Stationary increments**: The distribution of events depends only on the length of the interval. |
| 48 | +# 3. **Orderliness**: The probability of more than one event occurring in a small interval is negligible. |
| 49 | +# |
| 50 | +# Despite its simplicity, the homogeneous Poisson process provides a useful baseline model |
| 51 | +# and a reference point for more expressive point process models. Let's start simulating data using PointProcesses.jl |
| 52 | +using PointProcesses |
| 53 | +using Distributions |
| 54 | +using HypothesisTests |
| 55 | +using StatsBase |
| 56 | +using Plots |
| 57 | + |
| 58 | +using StableRNGs # hide |
| 59 | +rng = StableRNG(67); # hide |
| 60 | + |
| 61 | +# Usually we observe some data, typically, these will be event times. Assume, we have seen the following event times, from a cultured retinal ganglion cell neuron with no input stimulus: |
| 62 | +waiting_times = rand(rng, InverseGaussian(0.2, 0.5), 50) # hide |
| 63 | +event_times = cumsum(waiting_times) # hide |
| 64 | + |
| 65 | +h = History(event_times, 0.0, maximum(event_times) + 1.0) |
| 66 | + |
| 67 | +# Let's visualize our event times as a raster plot: |
| 68 | +raster = plot(; |
| 69 | + xlabel="Time", |
| 70 | + ylabel="Events", |
| 71 | + title="Raster Plot of Event Times", |
| 72 | + ylim=(-0.1, 1), |
| 73 | + yticks=false, |
| 74 | +) |
| 75 | +for time in h.times |
| 76 | + vline!([time]; color=:blue, linewidth=2, label=false) |
| 77 | +end |
| 78 | + |
| 79 | +raster |
| 80 | + |
| 81 | +# Looking at this plot, it's difficult to tell whether the event times are consistent with a homogeneous Poisson process. |
| 82 | +# We can bin the data in small time windows to get a sense of the event rate over time. From the assumptions above, we expect the distribution of |
| 83 | +# counts to be roughly Poisson distributed with rate λ * Δt in each bin of width Δt. |
| 84 | +bin_width = 1.0 |
| 85 | +bins = collect(h.tmin:bin_width:h.tmax) # bin edges |
| 86 | +counts = fit(Histogram, h.times, bins).weights; # counts per bin |
| 87 | + |
| 88 | +# Plot counts over time (counts/bin_width is a crude rate estimate) |
| 89 | +bin_centers = (bins[1:(end - 1)] .+ bins[2:end]) ./ 2 |
| 90 | + |
| 91 | +p_counts = plot( |
| 92 | + bin_centers, |
| 93 | + counts; |
| 94 | + seriestype=:bar, |
| 95 | + xlabel="Time", |
| 96 | + ylabel="Count per bin", |
| 97 | + title="Binned event counts", |
| 98 | +) |
| 99 | + |
| 100 | +p_counts |
| 101 | + |
| 102 | +# This is still difficult to tell. To get a better sense of whether the data is consistent with a homogeneous Poisson process, we can consider an alternative view of the homogeneous |
| 103 | +# point process based on the distribution of its waiting times. We discuss this below. |
| 104 | + |
| 105 | +# ### Two Equivalent Views |
| 106 | +# |
| 107 | +# The homogeneous Poisson process admits two equivalent and complementary descriptions: |
| 108 | +# |
| 109 | +# **(1) Inter-event times** |
| 110 | +# The waiting times between successive events are independent and identically distributed: |
| 111 | +# |
| 112 | +# τₖ = tₖ - tₖ₋₁ ∼ Exponential(λ) |
| 113 | +# |
| 114 | +# This view is particularly convenient for simulation. |
| 115 | +# |
| 116 | +# **(2) Counting process** |
| 117 | +# The total number of events in a time window of length T is distributed as: |
| 118 | +# |
| 119 | +# N(T) ∼ Poisson(λ T) |
| 120 | +# |
| 121 | +# Conditioned on N(T) = n, the event times are uniformly distributed on [0, T]. |
| 122 | +# |
| 123 | +# These two perspectives define the same stochastic process and will both be useful |
| 124 | +# for understanding likelihood-based inference. |
| 125 | + |
| 126 | +# Lets now calculate the waiting times and see how "exponential" they look: |
| 127 | +waiting_times = diff(vcat(h.tmin, h.times)); # includes waiting time from tmin to first event |
| 128 | + |
| 129 | +#- Rate estimate λ̂ = n / T |
| 130 | +T = duration(h) # should be h.tmax - h.tmin |
| 131 | +n = length(h) # number of events |
| 132 | +λ_est = n / T |
| 133 | + |
| 134 | +p_waiting = histogram( |
| 135 | + waiting_times; |
| 136 | + bins=10, |
| 137 | + normalize=:pdf, |
| 138 | + xlabel="Waiting time", |
| 139 | + ylabel="Density", |
| 140 | + title="Histogram of Waiting Times", |
| 141 | + label="Empirical", |
| 142 | +) |
| 143 | + |
| 144 | +x = range(0; stop=maximum(waiting_times), length=400) |
| 145 | +plot!( |
| 146 | + p_waiting, |
| 147 | + x, |
| 148 | + pdf.(Exponential(1 / λ_est), x); |
| 149 | + label="Exponential(λ̂=$(round(λ_est, digits=2)))", |
| 150 | + linewidth=2, |
| 151 | + color=:red, |
| 152 | +) |
| 153 | + |
| 154 | +p_waiting |
| 155 | + |
| 156 | +# This could be maybe exponential, but there is some signs that it is not. For example, why do we not see any waiting times below approximately 0.1s? |
| 157 | +# Let's fit the actual model now using PointProcesses.jl and run some formal statistical tests. |
| 158 | + |
| 159 | +# ### Fitting a Homogeneous Poisson Process |
| 160 | +pp_model = fit(PoissonProcess{Float64,Dirac{Nothing}}, h) |
| 161 | + |
| 162 | +println("Estimated rate λ̂ = ", pp_model.λ) # hide |
| 163 | + |
| 164 | +# We can now use this fitted model to run some statistical tests to see whether the homogeneous Poisson process is a good fit to the data. First let's plot a qq-plot |
| 165 | +# of the waiting times against the expected exponential distribution. |
| 166 | +q = range(0.01, 0.99; length=200) |
| 167 | +emp_q = quantile(waiting_times, q) |
| 168 | +theo_q = quantile.(Exponential(1 / λ_est), q) # Exp with mean 1/λ_est |
| 169 | + |
| 170 | +pqq = plot( |
| 171 | + theo_q, |
| 172 | + emp_q; |
| 173 | + seriestype=:scatter, |
| 174 | + xlabel="Theoretical quantiles: Exp(mean=1/λ̂)", |
| 175 | + ylabel="Empirical waiting-time quantiles", |
| 176 | + title="QQ plot: waiting times vs fitted exponential", |
| 177 | + aspect_ratio=:equal, |
| 178 | + label=false, |
| 179 | + xlim=(0, 1), |
| 180 | + ylim=(0, 1), |
| 181 | +) |
| 182 | +plot!(pqq, theo_q, theo_q; label=false) # y=x line |
| 183 | +pqq |
| 184 | + |
| 185 | +# The QQ-plot shows some deviations from the expected exponential distribution, particularly at the lower quantiles. We can use time-rescaling |
| 186 | +# to formally test the goodness-of-fit of the homogeneous Poisson process model. The time-rescaling theorem states that if we transform the event times tᵢ using the integrated |
| 187 | +# intensity function, the transformed times should be distributed as a homogeneous Poisson process with unit rate. Luckily PointProcesses.jl provides a convenient function to do this transformation for us: |
| 188 | +transformed_times = time_change(h, pp_model) |
| 189 | +transformed_waiting_times = diff(vcat(transformed_times.tmin, transformed_times.times)) |
| 190 | + |
| 191 | +# Under the null hypothesis that the data is generated by the fitted homogeneous Poisson process, these transformed waiting times should be i.i.d. Exp(1) distributed. |
| 192 | + |
| 193 | +ks_test = ExactOneSampleKSTest(transformed_waiting_times, Exponential(1.0)) |
| 194 | +println("KS test p-value: ", pvalue(ks_test)) # hide |
| 195 | + |
| 196 | +# The KS test p-value is very small, indicating that we can reject the null hypothesis that the data was generated by a homogeneous Poisson process at conventional significance levels. |
| 197 | +# This suggests that the homogeneous Poisson process is not an adequate model for the observed event times. There are several possible reasons for this, including: |
| 198 | +# - The event rate may not be constant over time (inhomogeneous Poisson process) |
| 199 | +# - There may be history dependence (e.g., refractory periods, bursting) |
| 200 | +# But this is a topic for another tutorial! |
0 commit comments