From e98f3e6cc5b4111866c973f05d1a2d88238a0999 Mon Sep 17 00:00:00 2001 From: jkling Date: Thu, 28 May 2026 19:44:11 +0200 Subject: [PATCH 1/5] Univariate Hawkes process - Added mark distribution to the process --- src/univariate/hawkes/fit.jl | 103 ++++++++++++ src/univariate/hawkes/hawkes_process.jl | 200 ++++-------------------- src/univariate/hawkes/simulation.jl | 37 +++++ test/hawkes.jl | 16 +- 4 files changed, 182 insertions(+), 174 deletions(-) create mode 100644 src/univariate/hawkes/fit.jl create mode 100644 src/univariate/hawkes/simulation.jl diff --git a/src/univariate/hawkes/fit.jl b/src/univariate/hawkes/fit.jl new file mode 100644 index 0000000..06c5922 --- /dev/null +++ b/src/univariate/hawkes/fit.jl @@ -0,0 +1,103 @@ +""" + StatsAPI.fit(rng, ::Type{HawkesProcess{T}}, h::History; step_tol::Float64 = 1e-6, max_iter::Int = 1000) where {T<:Real} + +Expectation-Maximization algorithm from [Lewis2011](@cite). +The relevant calculations are in page 4, equations 6-13. + +Let (t₁ < ... < tₙ) be the event times over the interval [0, T). We use the immigrant-descendant representation, +where immigrants arrive at a constant base rate μ and each each arrival may generate descendants following the +activation function α exp(-ω(t - tᵢ)). + +The algorithm consists in the following steps: +1. Start with some initial guess for the parameters μ, ψ, and ω. ψ = α ω is the branching factor. +2. Calculate λ(tᵢ; μ, ψ, ω) (`lambda_ts` in the code) using the procedure in [Ozaki1979](@cite). +3. For each tᵢ and each j < i, calculate Dᵢⱼ = P(tᵢ is a descendant of tⱼ) as + + Dᵢⱼ = ψ ω exp(-ω(tᵢ - tⱼ)) / λ(tᵢ; μ, ψ, ω). + + Define D = ∑_{j < i} Dᵢⱼ (expected number of descendants) and div = ∑_{j < i} (tᵢ - tⱼ) Dᵢⱼ. +4. Update the parameters as + μ = (N - D) / T + ψ = D / N + ω = D / div +5. If convergence criterion is met, return updated parameters, otherwise, back to step 2. + +Notice that, in the implementation, the process is normalized so the average inter-event time is equal to 1 and, +therefore, the interval of the process is transformed from T to N. Also, in equation (8) in the paper, + +∑_{i=1:n} pᵢᵢ = ∑_{i=1:n} (1 - ∑_{j < i} Dᵢⱼ) = N - D. + +""" +function StatsAPI.fit( + ::Type{HawkesProcess{T,MD}}, + h::History; + step_tol::Float64=1e-6, + max_iter::Int=1000, + rng::AbstractRNG=default_rng(), +) where {T<:Real,MD<:PointProcessMarkDistribution} + n = nb_events(h) + n == 0 && return HawkesProcess(zero(T), zero(T), zero(T)) + + tmax = T(duration(h)) + # Normalize times so average inter-event time is 1 (T -> n) + norm_ts = T.(h.times .* (n / tmax)) + + # preallocate + A = zeros(T, n) # A[i] = sum_{j= step_tol) && (n_iters < max_iter) + # compute A, S, and λ + for i in 2:n + Δ = norm_ts[i] - norm_ts[i - 1] + e = exp(-ω * Δ) + Ai_1 = A[i - 1] + A[i] = e * (one(T) + Ai_1) + S[i] = e * (S[i - 1] + Δ * (one(T) + Ai_1)) + lambda_ts[i] = μ + (ψ * ω) * A[i] + end + + # E-step aggregates + D = zero(T) # expected number of descendants + div = zero(T) # ∑ (t_i - t_j) D_{ij} + for i in 2:n + w = (ψ * ω) / lambda_ts[i] # factor common to all j= max_iter && @warn "Maximum number of iterations reached without convergence." + + mark_dist = fit(MD, h) + # Unnormalize back to original time scale (T -> tmax): + # parameters in normalized space (') relate to original by μ0=μ'*(n/tmax), ω0=ω'*(n/tmax), α0=(ψ'ω')*(n/tmax) + return HawkesProcess(μ * (n / tmax), ψ * ω * (n / tmax), ω * (n / tmax), mark_dist) +end + diff --git a/src/univariate/hawkes/hawkes_process.jl b/src/univariate/hawkes/hawkes_process.jl index 57ffad6..981c2c5 100644 --- a/src/univariate/hawkes/hawkes_process.jl +++ b/src/univariate/hawkes/hawkes_process.jl @@ -1,7 +1,7 @@ """ HawkesProcess{T<:Real} -Univariate Hawkes process with exponential decay kernel. +Univariate Marked Hawkes process with exponential decay kernel and events independent of marks. A Hawkes process is a self-exciting point process where each event increases the probability of future events. The conditional intensity function is given by: @@ -22,136 +22,55 @@ Conditions: Following the notation from [Lewis2011](@cite). """ -struct HawkesProcess{T<:Real} <: AbstractUnivariateProcess +struct HawkesProcess{T<:Real,D<:PointProcessMarkDistribution} <: AbstractUnivariateProcess μ::T α::T ω::T + mark_dist::D - function HawkesProcess(μ::T1, α::T2, ω::T3) where {T1,T2,T3} + function HawkesProcess(μ::T1, α::T2, ω::T3, mark_dist::D) where {T1,T2,T3<:Real,D<:PointProcessMarkDistribution} any((μ, α, ω) .< 0) && throw(DomainError((μ, α, ω), "All parameters must be non-negative.")) (α > 0 && α >= ω) && throw(DomainError((α, ω), "Parameter ω must be strictly smaller than α")) T = promote_type(T1, T2, T3) (μ_T, α_T, ω_T) = convert.(T, (μ, α, ω)) - new{T}(μ_T, α_T, ω_T) + new{T,D}(μ_T, α_T, ω_T, mark_dist) end end -function simulate(rng::AbstractRNG, hp::HawkesProcess, tmin, tmax) - sim = simulate_poisson_times(rng, hp.μ, tmin, tmax) # Simulate Poisson process with base rate - sim_desc = generate_descendants(rng, sim, tmax, hp.α, hp.ω) # Recursively generates descendants from first events - append!(sim, sim_desc) - sort!(sim) - return History(; times=sim, tmin=tmin, tmax=tmax, check_args=false) -end - -""" - StatsAPI.fit(rng, ::Type{HawkesProcess{T}}, h::History; step_tol::Float64 = 1e-6, max_iter::Int = 1000) where {T<:Real} - -Expectation-Maximization algorithm from [Lewis2011](@cite). -The relevant calculations are in page 4, equations 6-13. - -Let (t₁ < ... < tₙ) be the event times over the interval [0, T). We use the immigrant-descendant representation, -where immigrants arrive at a constant base rate μ and each each arrival may generate descendants following the -activation function α exp(-ω(t - tᵢ)). +HawkesProcess(μ, α, ω) = HawkesProcess(μ, α, ω, NoMarks()) -The algorithm consists in the following steps: -1. Start with some initial guess for the parameters μ, ψ, and ω. ψ = α ω is the branching factor. -2. Calculate λ(tᵢ; μ, ψ, ω) (`lambda_ts` in the code) using the procedure in [Ozaki1979](@cite). -3. For each tᵢ and each j < i, calculate Dᵢⱼ = P(tᵢ is a descendant of tⱼ) as +Base.ndims(hp::MultivariateHawkesProcess) = length(hp.μ) - Dᵢⱼ = ψ ω exp(-ω(tᵢ - tⱼ)) / λ(tᵢ; μ, ψ, ω). - - Define D = ∑_{j < i} Dᵢⱼ (expected number of descendants) and div = ∑_{j < i} (tᵢ - tⱼ) Dᵢⱼ. -4. Update the parameters as - μ = (N - D) / T - ψ = D / N - ω = D / div -5. If convergence criterion is met, return updated parameters, otherwise, back to step 2. - -Notice that, in the implementation, the process is normalized so the average inter-event time is equal to 1 and, -therefore, the interval of the process is transformed from T to N. Also, in equation (8) in the paper, - -∑_{i=1:n} pᵢᵢ = ∑_{i=1:n} (1 - ∑_{j < i} Dᵢⱼ) = N - D. +function ground_intensity(hp::HawkesProcess, h::History, t) + activation = sum(exp.(hp.ω .* (@view h.times[1:(searchsortedfirst(h.times, t) - 1)]))) + return hp.μ + (hp.α * activation / exp(hp.ω * t)) +end -""" -function StatsAPI.fit( - ::Type{HawkesProcess{T}}, - h::History; - step_tol::Float64=1e-6, - max_iter::Int=1000, - rng::AbstractRNG=default_rng(), -) where {T<:Real} - n = nb_events(h) - n == 0 && return HawkesProcess(zero(T), zero(T), zero(T)) - - tmax = T(duration(h)) - # Normalize times so average inter-event time is 1 (T -> n) - norm_ts = T.(h.times .* (n / tmax)) - - # preallocate - A = zeros(T, n) # A[i] = sum_{j= step_tol) && (n_iters < max_iter) - # compute A, S, and λ - for i in 2:n - Δ = norm_ts[i] - norm_ts[i - 1] - e = exp(-ω * Δ) - Ai_1 = A[i - 1] - A[i] = e * (one(T) + Ai_1) - S[i] = e * (S[i - 1] + Δ * (one(T) + Ai_1)) - lambda_ts[i] = μ + (ψ * ω) * A[i] - end - - # E-step aggregates - D = zero(T) # expected number of descendants - div = zero(T) # ∑ (t_i - t_j) D_{ij} - for i in 2:n - w = (ψ * ω) / lambda_ts[i] # factor common to all j= max_iter && @warn "Maximum number of iterations reached without convergence." - - # Unnormalize back to original time scale (T -> tmax): - # parameters in normalized space (') relate to original by μ0=μ'*(n/tmax), ω0=ω'*(n/tmax), α0=(ψ'ω')*(n/tmax) - return HawkesProcess(μ * (n / tmax), ψ * ω * (n / tmax), ω * (n / tmax)) + integral *= hp.α / hp.ω + integral += hp.μ * (tmax - tmin) # Integral of base rate + return integral end -# Type parameter for `HawkesProcess` was NOT explicitly provided -function StatsAPI.fit(HP::Type{HawkesProcess}, h::History{H,M}; kwargs...) where {H<:Real,M} - T = promote_type(Float64, H) - return fit(HP{T}, h; kwargs...) +function DensityInterface.logdensityof(hp::HawkesProcess, h::History) + A = zeros(nb_events(h)) # Vector A in Ozaki (1979) + for i in 2:nb_events(h) + A[i] = exp(-hp.ω * (h.times[i] - h.times[i - 1])) * (1 + A[i - 1]) + end + return sum(log.(hp.μ .+ (hp.α .* A))) - # Value of intensity at each event + (hp.μ * duration(h)) - # Integral of base rate + ((hp.α / hp.ω) * sum(1 .- exp.(-hp.ω .* (duration(h) .- h.times)))) + # Integral of each kernel + sum(log.([densityof(hp.mark_dist, t, h, m) for (t, m) in zip(h.times, h.marks)])) # Lok likelihood of marks end function time_change(h::History{R,M}, hp::HawkesProcess) where {R<:Real,M} @@ -180,60 +99,3 @@ function time_change(h::History{R,M}, hp::HawkesProcess) where {R<:Real,M} return History(; times=times, marks=h.marks, tmin=zero(T), tmax=tmax, check_args=false) # A time re-scaled process starts at t=0 end - -function ground_intensity(hp::HawkesProcess, h::History, t) - activation = sum(exp.(hp.ω .* (@view h.times[1:(searchsortedfirst(h.times, t) - 1)]))) - return hp.μ + (hp.α * activation / exp(hp.ω * t)) -end - -function integrated_ground_intensity(hp::HawkesProcess{T}, h::History, tmin, tmax) where {T} - U = promote_type(T, typeof(tmin), typeof(tmax)) - times = event_times(h, h.tmin, tmax) - integral = zero(U) - for ti in times - # Integral of activation function. 'max(tmin - ti, 0)' corrects for events that occurred - # inside or outside the interval [tmin, tmax]. - integral += (exp(-hp.ω * max(tmin - ti, 0)) - exp(-hp.ω * (tmax - ti))) - end - integral *= hp.α / hp.ω - integral += hp.μ * (tmax - tmin) # Integral of base rate - return integral -end - -function DensityInterface.logdensityof(hp::HawkesProcess, h::History) - A = zeros(nb_events(h)) # Vector A in Ozaki (1979) - for i in 2:nb_events(h) - A[i] = exp(-hp.ω * (h.times[i] - h.times[i - 1])) * (1 + A[i - 1]) - end - return sum(log.(hp.μ .+ (hp.α .* A))) - # Value of intensity at each event - (hp.μ * duration(h)) - # Integral of base rate - ((hp.α / hp.ω) * sum(1 .- exp.(-hp.ω .* (duration(h) .- h.times)))) # Integral of each kernel -end - -#= -Internal function for simulating Hawkes processes -The first generation, gen_0, is the `immigrants`, which is a set of event times. -For each t_g ∈ gen_n, simulate an inhomogeneous Poisson process over the interval [t_g, T] -with intensity λ(t) = α exp(-ω(t - t_g)) with the inverse method. -gen_{n+1} is the set of all events simulated from all events in gen_n. -The algorithm stops when the simulation from one generation results in no further events. -=# -function generate_descendants( - rng::AbstractRNG, immigrants::Vector{T}, tmax, α, ω -) where {T<:Real} - descendants = T[] - next_gen = immigrants - while !isempty(next_gen) - # OPTIMIZE: Can this be improved by avoiding allocations of `curr_gen` and `next_gen`? Or does the compiler take care of that? - curr_gen = copy(next_gen) # The current generation from which we simulate the next one - next_gen = eltype(immigrants)[] # Gathers all the descendants from the current generation - for parent in curr_gen # Generate the descendants for each individual event with the inverse method - activation_integral = (α / ω) * (one(T) - exp(ω * (parent - tmax))) - sim_transf = simulate_poisson_times(rng, one(T), zero(T), activation_integral) - @. sim_transf = parent - (inv(ω) * log(one(T) - ((ω / α) * sim_transf))) # Inverse of integral of the activation function - append!(next_gen, sim_transf) - end - append!(descendants, next_gen) - end - return descendants -end diff --git a/src/univariate/hawkes/simulation.jl b/src/univariate/hawkes/simulation.jl new file mode 100644 index 0000000..f436e92 --- /dev/null +++ b/src/univariate/hawkes/simulation.jl @@ -0,0 +1,37 @@ +function simulate(rng::AbstractRNG, hp::HawkesProcess, tmin, tmax) + sim = simulate_poisson_times(rng, hp.μ, tmin, tmax) # Simulate Poisson process with base rate + sim_desc = generate_descendants(rng, sim, tmax, hp.α, hp.ω) # Recursively generates descendants from first events + append!(sim, sim_desc) + sort!(sim) + h_temp = History(sim, tmin, tmax) + marks = [sample_mark(hp.mark_dist, t, h_temp) for t in event_times(h_temp)] + return History(; times=sim, tmin=tmin, tmax=tmax, marks=marks) +end + +#= +Internal function for simulating Hawkes processes +The first generation, gen_0, is the `immigrants`, which is a set of event times. +For each t_g ∈ gen_n, simulate an inhomogeneous Poisson process over the interval [t_g, T] +with intensity λ(t) = α exp(-ω(t - t_g)) with the inverse method. +gen_{n+1} is the set of all events simulated from all events in gen_n. +The algorithm stops when the simulation from one generation results in no further events. +=# +function generate_descendants( + rng::AbstractRNG, immigrants::Vector{T}, tmax, α, ω +) where {T<:Real} + descendants = T[] + next_gen = immigrants + while !isempty(next_gen) + # OPTIMIZE: Can this be improved by avoiding allocations of `curr_gen` and `next_gen`? Or does the compiler take care of that? + curr_gen = copy(next_gen) # The current generation from which we simulate the next one + next_gen = eltype(immigrants)[] # Gathers all the descendants from the current generation + for parent in curr_gen # Generate the descendants for each individual event with the inverse method + activation_integral = (α / ω) * (one(T) - exp(ω * (parent - tmax))) + sim_transf = simulate_poisson_times(rng, one(T), zero(T), activation_integral) + @. sim_transf = parent - (inv(ω) * log(one(T) - ((ω / α) * sim_transf))) # Inverse of integral of the activation function + append!(next_gen, sim_transf) + end + append!(descendants, next_gen) + end + return descendants +end diff --git a/test/hawkes.jl b/test/hawkes.jl index 009d058..05f4e77 100644 --- a/test/hawkes.jl +++ b/test/hawkes.jl @@ -1,6 +1,6 @@ # Constructor -@test HawkesProcess(1, 1, 2) isa HawkesProcess{Int} -@test HawkesProcess(1, 1, 2.0) isa HawkesProcess{Float64} +@test HawkesProcess(1, 1, 2) isa HawkesProcess{Int,NoMarks} +@test HawkesProcess(1, 1, 2.0, Normal()) isa HawkesProcess{Float64,Normal{Float64}} @test_throws DomainError HawkesProcess(1, 1, 1) @test_throws DomainError HawkesProcess(-1, 1, 2) @@ -39,17 +39,23 @@ h_sim = simulate(hp, 0.0, 10.0) @test isa(h_sim, History{Float64,Nothing}) @test isa(simulate(hp, BigFloat(0), BigFloat(10)), History{BigFloat,Nothing}) +hp_normal = HawkesProcess(1, 1, 2.0, Normal()) +h_sim2 = simulate(hp_normal, 0.0, 10.0) +@test issorted(h_sim2.times) +@test isa(h_sim2, History{Float64,Float64}) +@test isa(simulate(hp_normal, BigFloat(0), BigFloat(10)), History{BigFloat,Float64}) + # Fit Random.seed!(123) params_true = (100.0, 100.0, 200.0) model = HawkesProcess(params_true...) h_sim = simulate(model, 0.0, 50.0) -model_est = fit(HawkesProcess, h_sim) +model_est = fit(HawkesProcess{Float64,NoMarks}, h_sim) params_est = (model_est.μ, model_est.α, model_est.ω) @test isa(model_est, HawkesProcess) @test all((params_true .* 0.9) .<= params_est .<= (params_true .* 1.1)) -@test isa(fit(HawkesProcess, h_big), HawkesProcess{BigFloat}) -@test isa(fit(HawkesProcess{Float32}, h_big), HawkesProcess{Float32}) +@test isa(fit(HawkesProcess{BigFloat,NoMarks}, h_big), HawkesProcess{BigFloat}) +@test isa(fit(HawkesProcess{Float64,Normal}, h_sim2), HawkesProcess{Float64,Normal{Float64}}) # logdensityof @test logdensityof(hp, h) ≈ From 7f45b10ed1d8f6c59c2982a1efa413206844b9a0 Mon Sep 17 00:00:00 2001 From: jkling Date: Tue, 2 Jun 2026 22:43:14 +0200 Subject: [PATCH 2/5] Interface for Multivariate PointProcesses - Added `univariate_process.jl` and `multivariate_process.jl` files defining the common interfaces - Updated the multivariate processes to comply with the interface - Split univariate Hawkes implementation into thre files --- src/HypothesisTests/PPTests/bootstrap_test.jl | 2 +- .../PPTests/monte_carlo_test.jl | 2 +- src/HypothesisTests/point_process_tests.jl | 4 +- src/PointProcesses.jl | 4 ++ src/abstract_point_process.jl | 45 +++++-------------- src/history.jl | 5 --- src/multivariate/independent_multivariate.jl | 36 +++------------ src/multivariate_process.jl | 41 +++++++++++++++++ src/univariate/hawkes/hawkes_process.jl | 2 - src/univariate_process.jl | 31 +++++++++++++ test/hawkes.jl | 3 +- 11 files changed, 97 insertions(+), 78 deletions(-) create mode 100644 src/multivariate_process.jl create mode 100644 src/univariate_process.jl diff --git a/src/HypothesisTests/PPTests/bootstrap_test.jl b/src/HypothesisTests/PPTests/bootstrap_test.jl index 228f65f..1eacdbf 100644 --- a/src/HypothesisTests/PPTests/bootstrap_test.jl +++ b/src/HypothesisTests/PPTests/bootstrap_test.jl @@ -71,7 +71,7 @@ function BootstrapTest( end # Estimate process and calculate statistic from data - pp_est = fit(PP, h; rng=rng) + pp_est = fit(PP, h; rng=rng)::AbstractPointProcess # JET complains if not specified stat = statistic(S, pp_est, h) # Initialize vector with test statistics from simulations diff --git a/src/HypothesisTests/PPTests/monte_carlo_test.jl b/src/HypothesisTests/PPTests/monte_carlo_test.jl index d1fc8d7..26e27fb 100644 --- a/src/HypothesisTests/PPTests/monte_carlo_test.jl +++ b/src/HypothesisTests/PPTests/monte_carlo_test.jl @@ -110,6 +110,6 @@ function MonteCarloTest( throw(ArgumentError("Test is not valid for empty event history.")) end - pp_est = fit(PP, h) + pp_est = fit(PP, h)::AbstractPointProcess # JET complains if not specified return MonteCarloTest(S, pp_est, h; kwargs...) end diff --git a/src/HypothesisTests/point_process_tests.jl b/src/HypothesisTests/point_process_tests.jl index 1f4d0eb..d268e20 100644 --- a/src/HypothesisTests/point_process_tests.jl +++ b/src/HypothesisTests/point_process_tests.jl @@ -26,12 +26,12 @@ end Internal function for performing the appropriate transformation on the event times according to the selected distribution. =# -function transform(::Type{<:Uniform}, pp::AbstractPointProcess, h) +function transform(::Type{<:Uniform}, pp::AbstractPointProcess, h::History) transf = time_change(h, pp) # transf → time re-scaled event times return transf.times, Uniform(transf.tmin, transf.tmax) end -function transform(::Type{<:Exponential}, pp::AbstractPointProcess, h) +function transform(::Type{<:Exponential}, pp::AbstractPointProcess, h::History) inter_transf = diff(time_change(h, pp).times) # inter_transf → sorted time transformed inter event times sort!(inter_transf) return inter_transf, Exponential(1) diff --git a/src/PointProcesses.jl b/src/PointProcesses.jl index 5a7d717..81f97fb 100644 --- a/src/PointProcesses.jl +++ b/src/PointProcesses.jl @@ -87,6 +87,8 @@ export PointProcessTest, BootstrapTest, MonteCarloTest include("history.jl") include("mark_distributions.jl") include("abstract_point_process.jl") +include("univariate_process.jl") +include("multivariate_process.jl") include("simulation.jl") include("bounded_point_process.jl") @@ -104,6 +106,8 @@ include("univariate/poisson/inhomogeneous/fit.jl") ### Hawkes include("univariate/hawkes/hawkes_process.jl") +include("univariate/hawkes/fit.jl") +include("univariate/hawkes/simulation.jl") ## Multivariate processes diff --git a/src/abstract_point_process.jl b/src/abstract_point_process.jl index 263d154..2fe152f 100644 --- a/src/abstract_point_process.jl +++ b/src/abstract_point_process.jl @@ -5,20 +5,6 @@ Common interface for all temporal point processes. """ abstract type AbstractPointProcess end -""" - AbstractUnivariateProcess - -Abstract type for univariate temporal point processes. -""" -abstract type AbstractUnivariateProcess <: AbstractPointProcess end - -""" - AbstractMultivariateProcess - -Abstract type for multivariate temporal point processes. -""" -abstract type AbstractMultivariateProcess <: AbstractPointProcess end - @inline DensityInterface.DensityKind(::AbstractPointProcess) = HasDensity() """ @@ -28,10 +14,7 @@ Return the number of dimensions for a temporal point process `pp`. """ Base.ndims(::AbstractPointProcess) -Base.ndims(::AbstractUnivariateProcess) = 1 - ## Intensity functions - """ ground_intensity(pp, h, t) @@ -55,7 +38,7 @@ For multivariate processes, it returns a vector of mark distributions for each d mark_distribution(pp, h, t, d) computes the mark distribution for a multivariate process `pp` at dimension `d`. """ -mark_distribution(pp::AbstractPointProcess, t, h) = mark_distribution(pp.mark_dist, t, h) +function mark_distribution end """ intensity(pp, m, t, h) @@ -67,9 +50,7 @@ intensity(pp, h, t, d) computes the intensity for a multivariate process `pp` at The conditional intensity function `λ(t,m|h)` quantifies the instantaneous risk of an event with mark `m` occurring at time `t` after history `h`. """ -function intensity(pp::AbstractPointProcess, m, t, h) - return ground_intensity(pp, t, h) * densityof(pp.mark_dist, t, h, m) -end +function intensity end """ log_intensity(pp, m, t, h) @@ -79,9 +60,7 @@ For multivariate processes, it returns a vector of log intensities for each dime log_intensity(pp, h, t, d) computes the log intensity for a multivariate process `pp` at dimension `d`. """ -function log_intensity(pp::AbstractPointProcess, m, t, h) - return log(intensity(pp, m, t, h)) -end +function log_intensity end ## Simulation @@ -122,13 +101,7 @@ Compute the log probability density function for a temporal point process `pp` a ``` The default method uses a loop over events combined with `integrated_ground_intensity`, but it should be reimplemented for specific processes if faster computation is possible. """ -function DensityInterface.logdensityof(pp::AbstractPointProcess, h::History) - l = -integrated_ground_intensity(pp, h, min_time(h), max_time(h)) - for (t, m) in zip(event_times(h), event_marks(h)) - l += log_intensity(pp, m, t, h) - end - return l -end +DensityInterface.logdensityof """ fit(::Type{PP}, h) @@ -150,7 +123,9 @@ Not implemented by default. """ function fit_map end -function time_change(h::History, pp::AbstractPointProcess) - Λ(t) = integrated_ground_intensity(pp, h, min_time(h), t) - return time_change(h, Λ) -end +""" + time_change + +Apply the time rescaling `t -> Λ(t)` to history `h`. +""" +function time_change end \ No newline at end of file diff --git a/src/history.jl b/src/history.jl index be761f0..17eb91e 100644 --- a/src/history.jl +++ b/src/history.jl @@ -444,11 +444,6 @@ function Base.cat(h1::History, h2::History) ) end -""" - time_change(h, Λ) - -Apply the time rescaling `t -> Λ(t)` to history `h`. -""" function time_change(h::History, Λ) new_times = Λ.(event_times(h)) new_marks = copy(event_marks(h)) diff --git a/src/multivariate/independent_multivariate.jl b/src/multivariate/independent_multivariate.jl index ff680c3..b9f21a4 100644 --- a/src/multivariate/independent_multivariate.jl +++ b/src/multivariate/independent_multivariate.jl @@ -18,56 +18,30 @@ end Base.ndims(pp::IndependentMultivariateProcess) = length(pp.processes) ## AbstractPointProcess interface -function ground_intensity(pp::IndependentMultivariateProcess, t, h, d) +function ground_intensity(pp::IndependentMultivariateProcess, t, h::History, d) return ground_intensity(pp.processes[d], t, History(h, d)) end -function ground_intensity(pp::IndependentMultivariateProcess, t, h) - return [ground_intensity(pp, t, h, d) for d in 1:ndims(pp)] +function intensity(pp::IndependentMultivariateProcess, m, t, h::History, d) + return intensity(pp.processes[d], m, t, History(h, d)) end -function mark_distribution(pp::IndependentMultivariateProcess, t, h, d) +function mark_distribution(pp::IndependentMultivariateProcess, t, h::History, d) return mark_distribution(pp.processes[d], t, History(h, d)) end -function mark_distribution(pp::IndependentMultivariateProcess, t, h) - return [mark_distribution(pp, t, h, d) for d in 1:ndims(pp)] -end - function mark_distribution(pp::IndependentMultivariateProcess, t) return [mark_distribution(pp.processes[d], t) for d in 1:ndims(pp)] end -function intensity(pp::IndependentMultivariateProcess, m, t, h, d) - return intensity(pp.processes[d], m, t, History(h, d)) -end - -function intensity(pp::IndependentMultivariateProcess, m, t, h) - return [intensity(pp, m, t, h, d) for d in 1:ndims(pp)] -end - -function log_intensity(pp::IndependentMultivariateProcess, m, t, h, d) - log(intensity(pp, m, t, h, d)) -end - -log_intensity(pp::IndependentMultivariateProcess, m, t, h) = log.(intensity(pp, m, t, h)) - -function ground_intensity_bound(pp::IndependentMultivariateProcess, t, h, d) +function ground_intensity_bound(pp::IndependentMultivariateProcess, t, h::History, d) return ground_intensity_bound(pp.processes[d], t, History(h, d)) end -function ground_intensity_bound(pp::IndependentMultivariateProcess, t, h) - return [ground_intensity_bound(pp, t, h, d) for d in 1:ndims(pp)] -end - function integrated_ground_intensity(pp::IndependentMultivariateProcess, h, a, b, d) return integrated_ground_intensity(pp.processes[d], History(h, d), a, b) end -function integrated_ground_intensity(pp::IndependentMultivariateProcess, h, a, b) - return [integrated_ground_intensity(pp, h, a, b, d) for d in 1:ndims(pp)] -end - function DensityInterface.logdensityof(pp::IndependentMultivariateProcess, h::History) return sum( logdensityof( diff --git a/src/multivariate_process.jl b/src/multivariate_process.jl new file mode 100644 index 0000000..1b7f40e --- /dev/null +++ b/src/multivariate_process.jl @@ -0,0 +1,41 @@ +""" + AbstractMultivariateProcess + +Abstract type for multivariate temporal point processes. +""" +abstract type AbstractMultivariateProcess <: AbstractPointProcess end + +Base.ndims(pp::AbstractMultivariateProcess) = length(pp.mark_dist) + +function mark_distribution(pp::AbstractMultivariateProcess, t, h) + return [mark_distribution(pp, t, h, d) for d in 1:ndims(pp)] +end + +function ground_intensity(pp::AbstractMultivariateProcess, t, h::History) + return [ground_intensity(pp, t, h, d) for d in 1:ndims(pp)] +end + +function integrated_ground_intensity(pp::AbstractMultivariateProcess, h::History, a, b) + return [integrated_ground_intensity(pp, h, a, b, d) for d in 1:ndims(pp)] +end + +function ground_intensity_bound(pp::AbstractMultivariateProcess, t, h::History) + return [ground_intensity_bound(pp, t, h, d) for d in 1:ndims(pp)] +end + +function intensity(pp::AbstractMultivariateProcess, m, t, h, d) + return ground_intensity(pp, t, h, d) * densityof(pp.mark_dist[d], t, h, m) +end + +function intensity(pp::AbstractMultivariateProcess, m, t, h::History) + return [intensity(pp, m, t, h, d) for d in 1:ndims(pp)] +end + +function log_intensity(pp::AbstractMultivariateProcess, m, t, h, d) + return log(intensity(pp, m, t, h, d)) +end + +function log_intensity(pp::AbstractMultivariateProcess, m, t, h::History) + return [log_intensity(pp, m, t, h, d) for d in 1:ndims(pp)] +end + diff --git a/src/univariate/hawkes/hawkes_process.jl b/src/univariate/hawkes/hawkes_process.jl index 981c2c5..9cfe382 100644 --- a/src/univariate/hawkes/hawkes_process.jl +++ b/src/univariate/hawkes/hawkes_process.jl @@ -41,8 +41,6 @@ end HawkesProcess(μ, α, ω) = HawkesProcess(μ, α, ω, NoMarks()) -Base.ndims(hp::MultivariateHawkesProcess) = length(hp.μ) - function ground_intensity(hp::HawkesProcess, h::History, t) activation = sum(exp.(hp.ω .* (@view h.times[1:(searchsortedfirst(h.times, t) - 1)]))) return hp.μ + (hp.α * activation / exp(hp.ω * t)) diff --git a/src/univariate_process.jl b/src/univariate_process.jl new file mode 100644 index 0000000..eb15ded --- /dev/null +++ b/src/univariate_process.jl @@ -0,0 +1,31 @@ +""" + AbstractUnivariateProcess + +Abstract type for univariate temporal point processes. +""" +abstract type AbstractUnivariateProcess <: AbstractPointProcess end + +Base.ndims(::AbstractUnivariateProcess) = 1 + +mark_distribution(pp::AbstractUnivariateProcess, t, h) = mark_distribution(pp.mark_dist, t, h) + +function intensity(pp::AbstractUnivariateProcess, m, t, h) + return ground_intensity(pp, t, h) * densityof(pp.mark_dist, t, h, m) +end + +function log_intensity(pp::AbstractUnivariateProcess, m, t, h) + return log(intensity(pp, m, t, h)) +end + +function DensityInterface.logdensityof(pp::AbstractUnivariateProcess, h::History) + l = -integrated_ground_intensity(pp, h, min_time(h), max_time(h)) + for (t, m) in zip(event_times(h), event_marks(h)) + l += log_intensity(pp, m, t, h) + end + return l +end + +function time_change(h::History, pp::AbstractUnivariateProcess) + Λ(t) = integrated_ground_intensity(pp, h, min_time(h), t) + return time_change(h, Λ) +end \ No newline at end of file diff --git a/test/hawkes.jl b/test/hawkes.jl index 05f4e77..0e5eb0b 100644 --- a/test/hawkes.jl +++ b/test/hawkes.jl @@ -58,7 +58,8 @@ params_est = (model_est.μ, model_est.α, model_est.ω) @test isa(fit(HawkesProcess{Float64,Normal}, h_sim2), HawkesProcess{Float64,Normal{Float64}}) # logdensityof -@test logdensityof(hp, h) ≈ +h_nomarks = History(h.times, h.tmin, h.tmax) +@test logdensityof(hp, h_nomarks) ≈ sum(log.(hp.μ .+ (hp.α .* [0, exp(-hp.ω), exp(-hp.ω * 2) + exp(-hp.ω * 3)]))) - integral From c2d4ba1c356aaaa6ea2324e4c714758b19730215 Mon Sep 17 00:00:00 2001 From: jkling Date: Sat, 20 Jun 2026 10:59:54 +0200 Subject: [PATCH 3/5] Automatic formatting --- src/multivariate_process.jl | 1 - src/univariate/hawkes/fit.jl | 1 - src/univariate/hawkes/hawkes_process.jl | 4 +++- src/univariate_process.jl | 6 ++++-- test/hawkes.jl | 4 +++- 5 files changed, 10 insertions(+), 6 deletions(-) diff --git a/src/multivariate_process.jl b/src/multivariate_process.jl index 1b7f40e..67e662f 100644 --- a/src/multivariate_process.jl +++ b/src/multivariate_process.jl @@ -38,4 +38,3 @@ end function log_intensity(pp::AbstractMultivariateProcess, m, t, h::History) return [log_intensity(pp, m, t, h, d) for d in 1:ndims(pp)] end - diff --git a/src/univariate/hawkes/fit.jl b/src/univariate/hawkes/fit.jl index 06c5922..c4ec8b6 100644 --- a/src/univariate/hawkes/fit.jl +++ b/src/univariate/hawkes/fit.jl @@ -100,4 +100,3 @@ function StatsAPI.fit( # parameters in normalized space (') relate to original by μ0=μ'*(n/tmax), ω0=ω'*(n/tmax), α0=(ψ'ω')*(n/tmax) return HawkesProcess(μ * (n / tmax), ψ * ω * (n / tmax), ω * (n / tmax), mark_dist) end - diff --git a/src/univariate/hawkes/hawkes_process.jl b/src/univariate/hawkes/hawkes_process.jl index 9cfe382..52079e8 100644 --- a/src/univariate/hawkes/hawkes_process.jl +++ b/src/univariate/hawkes/hawkes_process.jl @@ -28,7 +28,9 @@ struct HawkesProcess{T<:Real,D<:PointProcessMarkDistribution} <: AbstractUnivari ω::T mark_dist::D - function HawkesProcess(μ::T1, α::T2, ω::T3, mark_dist::D) where {T1,T2,T3<:Real,D<:PointProcessMarkDistribution} + function HawkesProcess( + μ::T1, α::T2, ω::T3, mark_dist::D + ) where {T1,T2,T3<:Real,D<:PointProcessMarkDistribution} any((μ, α, ω) .< 0) && throw(DomainError((μ, α, ω), "All parameters must be non-negative.")) (α > 0 && α >= ω) && diff --git a/src/univariate_process.jl b/src/univariate_process.jl index eb15ded..4ad7c44 100644 --- a/src/univariate_process.jl +++ b/src/univariate_process.jl @@ -7,7 +7,9 @@ abstract type AbstractUnivariateProcess <: AbstractPointProcess end Base.ndims(::AbstractUnivariateProcess) = 1 -mark_distribution(pp::AbstractUnivariateProcess, t, h) = mark_distribution(pp.mark_dist, t, h) +function mark_distribution(pp::AbstractUnivariateProcess, t, h) + mark_distribution(pp.mark_dist, t, h) +end function intensity(pp::AbstractUnivariateProcess, m, t, h) return ground_intensity(pp, t, h) * densityof(pp.mark_dist, t, h, m) @@ -28,4 +30,4 @@ end function time_change(h::History, pp::AbstractUnivariateProcess) Λ(t) = integrated_ground_intensity(pp, h, min_time(h), t) return time_change(h, Λ) -end \ No newline at end of file +end diff --git a/test/hawkes.jl b/test/hawkes.jl index 0e5eb0b..65a277a 100644 --- a/test/hawkes.jl +++ b/test/hawkes.jl @@ -55,7 +55,9 @@ params_est = (model_est.μ, model_est.α, model_est.ω) @test isa(model_est, HawkesProcess) @test all((params_true .* 0.9) .<= params_est .<= (params_true .* 1.1)) @test isa(fit(HawkesProcess{BigFloat,NoMarks}, h_big), HawkesProcess{BigFloat}) -@test isa(fit(HawkesProcess{Float64,Normal}, h_sim2), HawkesProcess{Float64,Normal{Float64}}) +@test isa( + fit(HawkesProcess{Float64,Normal}, h_sim2), HawkesProcess{Float64,Normal{Float64}} +) # logdensityof h_nomarks = History(h.times, h.tmin, h.tmax) From af47ed79e26066ccba6f30e3ca41258deb0b2f25 Mon Sep 17 00:00:00 2001 From: jkling Date: Sun, 21 Jun 2026 18:28:49 +0200 Subject: [PATCH 4/5] Fixed convergence in example `CustomProcesses` --- docs/examples/CustomProcesses.jl | 11 +++++++++-- docs/examples/Hawkes.jl | 2 +- src/univariate/hawkes/hawkes_process.jl | 2 +- src/univariate_process.jl | 2 +- 4 files changed, 12 insertions(+), 5 deletions(-) diff --git a/docs/examples/CustomProcesses.jl b/docs/examples/CustomProcesses.jl index c8119e5..2d0c9c2 100644 --- a/docs/examples/CustomProcesses.jl +++ b/docs/examples/CustomProcesses.jl @@ -34,7 +34,6 @@ using Distributions using Plots using StatsAPI using Optim - import PointProcesses: NoMarks, AbstractMarkDistribution, AbstractUnivariateProcess struct TwoStateModel{R<:Real,D<:PointProcessMarkDistribution} <: AbstractUnivariateProcess @@ -147,7 +146,15 @@ function StatsAPI.fit( lower_bound = [0.0, 0.0, 0.0] upper_bound = [Inf, Inf, Inf] - result = optimize(objective, lower_bound, upper_bound, init_params) + result = optimize( + objective, + lower_bound, + upper_bound, + init_params, + NelderMead(), + Optim.Options(; x_reltol=1e-3), + ) + optimal = Optim.minimizer(result) return TwoStateModel(optimal..., NoMarks()) diff --git a/docs/examples/Hawkes.jl b/docs/examples/Hawkes.jl index 59a810f..80d887b 100644 --- a/docs/examples/Hawkes.jl +++ b/docs/examples/Hawkes.jl @@ -126,7 +126,7 @@ plot( # The goal of this analysis is to understand the self-exciting nature of the litter box entries. I.e., if one cat uses the litter box, # does that increase the likelihood of another cat using it soon after? To do this, we can use the implementation in PointProcesses.jl full_history = History(sort(data.t), 0.0, maximum(data.t) + 1.0) -hawkes_model = fit(HawkesProcess, full_history) +hawkes_model = fit(HawkesProcess{Float64,NoMarks}, full_history) println("Fitted Hawkes Process Parameters:") # hide println("Base intensity (μ): ", hawkes_model.μ) # hide diff --git a/src/univariate/hawkes/hawkes_process.jl b/src/univariate/hawkes/hawkes_process.jl index 52079e8..f806862 100644 --- a/src/univariate/hawkes/hawkes_process.jl +++ b/src/univariate/hawkes/hawkes_process.jl @@ -37,7 +37,7 @@ struct HawkesProcess{T<:Real,D<:PointProcessMarkDistribution} <: AbstractUnivari throw(DomainError((α, ω), "Parameter ω must be strictly smaller than α")) T = promote_type(T1, T2, T3) (μ_T, α_T, ω_T) = convert.(T, (μ, α, ω)) - new{T,D}(μ_T, α_T, ω_T, mark_dist) + return new{T,D}(μ_T, α_T, ω_T, mark_dist) end end diff --git a/src/univariate_process.jl b/src/univariate_process.jl index 4ad7c44..fd91d85 100644 --- a/src/univariate_process.jl +++ b/src/univariate_process.jl @@ -8,7 +8,7 @@ abstract type AbstractUnivariateProcess <: AbstractPointProcess end Base.ndims(::AbstractUnivariateProcess) = 1 function mark_distribution(pp::AbstractUnivariateProcess, t, h) - mark_distribution(pp.mark_dist, t, h) + return mark_distribution(pp.mark_dist, t, h) end function intensity(pp::AbstractUnivariateProcess, m, t, h) From 4052588218ec88afb7afda45fec057a752600af7 Mon Sep 17 00:00:00 2001 From: jkling Date: Thu, 2 Jul 2026 17:57:01 +0200 Subject: [PATCH 5/5] PR fixes Added tests for custom custom processes --- docs/examples/CustomProcesses.jl | 52 ++++++++++---- src/HypothesisTests/PPTests/bootstrap_test.jl | 2 +- src/abstract_point_process.jl | 16 ----- src/bounded_point_process.jl | 9 +++ src/history.jl | 11 ++- src/multivariate/independent_multivariate.jl | 16 ++--- src/multivariate_process.jl | 36 +++++++++- src/simulation.jl | 4 ++ src/univariate/hawkes/fit.jl | 10 +-- src/univariate/hawkes/hawkes_process.jl | 14 ++-- src/univariate/hawkes/simulation.jl | 8 ++- src/univariate/poisson/fit.jl | 11 ++- .../inhomogeneous_poisson_process.jl | 4 +- src/univariate_process.jl | 22 ++++-- test/bounded_point_process.jl | 17 ----- test/custom_processes.jl | 72 +++++++++++++++++++ test/hawkes.jl | 4 +- test/runtests.jl | 3 + 18 files changed, 222 insertions(+), 89 deletions(-) create mode 100644 test/custom_processes.jl diff --git a/docs/examples/CustomProcesses.jl b/docs/examples/CustomProcesses.jl index 2d0c9c2..a7a96c1 100644 --- a/docs/examples/CustomProcesses.jl +++ b/docs/examples/CustomProcesses.jl @@ -33,6 +33,7 @@ using PointProcesses using Distributions using Plots using StatsAPI +using StatsBase using Optim import PointProcesses: NoMarks, AbstractMarkDistribution, AbstractUnivariateProcess @@ -140,23 +141,45 @@ plot( # process parameters. function StatsAPI.fit( - ::Type{TwoStateModel{R1,NoMarks}}, h::History, init_params::Vector{R2} -) where {R1<:Real,R2<:Real} + ::Type{TwoStateModel{R1,NoMarks}}, h::History; max_tries::Int=100 +) where {R1<:Real} objective(params) = -logdensityof(TwoStateModel(params..., NoMarks()), h) lower_bound = [0.0, 0.0, 0.0] upper_bound = [Inf, Inf, Inf] - result = optimize( - objective, - lower_bound, - upper_bound, - init_params, - NelderMead(), - Optim.Options(; x_reltol=1e-3), - ) - - optimal = Optim.minimizer(result) - + mean_intensity = nb_events(h) / duration(h) + mean_interarrival = mean(diff(h.times)) + solved = false + n_tries = 1 + optimal = nothing + + while !solved && n_tries <= max_tries + try # `Optmi.jl` may call the objective function with parameters outside of the constraints + init_params = [ + 0.25 * mean_intensity + rand() * 0.5 * mean_intensity, + 0.25 * mean_intensity + rand() * 0.5 * mean_intensity, + 0.1 * mean_interarrival + rand() * 0.8 * mean_interarrival, + ] + + result = optimize( + objective, + lower_bound, + upper_bound, + init_params, + NelderMead(), + Optim.Options(; x_reltol=1e-3), + ) + + optimal = Optim.minimizer(result) + solved = true + + catch DomainError + n_tries += 1 + end + end + if !solved + throw(ErrorException("Solver did not converge.")) + end return TwoStateModel(optimal..., NoMarks()) end @@ -175,8 +198,7 @@ end true_params = random_params() h_unknown = simulate(TwoStateModel(true_params..., NoMarks()), 0.0, 100.0) -init_params = random_params() -pp_estimated = fit(TwoStateModel{Float64,NoMarks}, h_unknown, init_params) +pp_estimated = fit(TwoStateModel{Float64,NoMarks}, h_unknown) estimated_params = [pp_estimated.λ, pp_estimated.Δ, pp_estimated.τ] # hide println("True and estimated parameters:") # hide diff --git a/src/HypothesisTests/PPTests/bootstrap_test.jl b/src/HypothesisTests/PPTests/bootstrap_test.jl index 12b8346..839d5a3 100644 --- a/src/HypothesisTests/PPTests/bootstrap_test.jl +++ b/src/HypothesisTests/PPTests/bootstrap_test.jl @@ -71,7 +71,7 @@ function BootstrapTest( end # Estimate process and calculate statistic from data - pp_est = fit(PP, h; rng=rng)::AbstractPointProcess # JET complains if not specified + pp_est = fit(PP, h; rng=rng) stat = statistic(S, pp_est, h) # Initialize vector with test statistics from simulations diff --git a/src/abstract_point_process.jl b/src/abstract_point_process.jl index 77e8141..8003676 100644 --- a/src/abstract_point_process.jl +++ b/src/abstract_point_process.jl @@ -92,17 +92,6 @@ integrated_ground_intensity(pp, h, a, b, d) computes the integrated ground inten """ function integrated_ground_intensity end -""" - logdensityof(pp, h) - -Compute the log probability density function for a temporal point process `pp` applied to history `h`: -``` -ℓ(h) = Σₖ log λ(tₖ|hₖ) - Λ(h) -``` -The default method uses a loop over events combined with `integrated_ground_intensity`, but it should be reimplemented for specific processes if faster computation is possible. -""" -DensityInterface.logdensityof - """ fit(::Type{PP}, h) fit(::Type{PP}, histories) @@ -123,11 +112,6 @@ Not implemented by default. """ function fit_map end -function time_change(h::History, pp::AbstractPointProcess) - Λ(t) = integrated_ground_intensity(pp, h, min_time(h), t) - return time_change(h, Λ) -end - """ simulate([rng,] pp, tmin, tmax) diff --git a/src/bounded_point_process.jl b/src/bounded_point_process.jl index 057832a..c123318 100644 --- a/src/bounded_point_process.jl +++ b/src/bounded_point_process.jl @@ -39,3 +39,12 @@ end function integrated_ground_intensity(bpp::BoundedPointProcess, args...) return integrated_ground_intensity(bpp.pp, args...) end + +function time_change(h::History{T}, pp::BoundedPointProcess) where {T} + if min_time(h) < min_time(pp) || max_time(h) > max_time(pp) + throw(ArgumentError("History is defined outside the bounds of the process")) + end + return time_change(h, pp.pp) +end + +simulate(pp::BoundedPointProcess) = simulate(pp.pp, pp.tmin, pp.tmax) diff --git a/src/history.jl b/src/history.jl index 79b916f..d0ef40e 100644 --- a/src/history.jl +++ b/src/history.jl @@ -153,8 +153,13 @@ end function History(; times, tmin, tmax, marks=nothing, dims=nothing, check_args=true) if times isa Vector{<:Real} marks === nothing && (marks = fill(nothing, length(times))) - dims === nothing && (dims = fill(nothing, length(times))) - return History(times, tmin, tmax, marks, dims; check_args=check_args) + if dims === nothing + dims = fill(nothing, length(times)) + N = 1 + else + N = length(unique(dims)) + end + return History(times, tmin, tmax, marks, dims, N; check_args=check_args) else marks === nothing && (marks = [fill(nothing, length(times[i])) for i in 1:length(times)]) @@ -464,7 +469,7 @@ function Base.cat(h1::History, h2::History) ) end -function time_change(h::History, Λ) +function time_change(h::History{T}, Λ) where {T} new_times = Λ.(event_times(h)) new_marks = copy(event_marks(h)) new_tmin = Λ(min_time(h)) diff --git a/src/multivariate/independent_multivariate.jl b/src/multivariate/independent_multivariate.jl index b9f21a4..aff4ddd 100644 --- a/src/multivariate/independent_multivariate.jl +++ b/src/multivariate/independent_multivariate.jl @@ -18,27 +18,25 @@ end Base.ndims(pp::IndependentMultivariateProcess) = length(pp.processes) ## AbstractPointProcess interface -function ground_intensity(pp::IndependentMultivariateProcess, t, h::History, d) +function ground_intensity(pp::IndependentMultivariateProcess, t, h::History, d::Int) return ground_intensity(pp.processes[d], t, History(h, d)) end -function intensity(pp::IndependentMultivariateProcess, m, t, h::History, d) +function intensity(pp::IndependentMultivariateProcess, m, t, h::History, d::Int) return intensity(pp.processes[d], m, t, History(h, d)) end -function mark_distribution(pp::IndependentMultivariateProcess, t, h::History, d) +function mark_distribution(pp::IndependentMultivariateProcess, t, h::History, d::Int) return mark_distribution(pp.processes[d], t, History(h, d)) end -function mark_distribution(pp::IndependentMultivariateProcess, t) - return [mark_distribution(pp.processes[d], t) for d in 1:ndims(pp)] -end - -function ground_intensity_bound(pp::IndependentMultivariateProcess, t, h::History, d) +function ground_intensity_bound(pp::IndependentMultivariateProcess, t, h::History, d::Int) return ground_intensity_bound(pp.processes[d], t, History(h, d)) end -function integrated_ground_intensity(pp::IndependentMultivariateProcess, h, a, b, d) +function integrated_ground_intensity( + pp::IndependentMultivariateProcess, h::History, a, b, d::Int +) return integrated_ground_intensity(pp.processes[d], History(h, d), a, b) end diff --git a/src/multivariate_process.jl b/src/multivariate_process.jl index 67e662f..a6425b7 100644 --- a/src/multivariate_process.jl +++ b/src/multivariate_process.jl @@ -2,12 +2,26 @@ AbstractMultivariateProcess Abstract type for multivariate temporal point processes. + +To implement a multivariate process, one must subtype `AbstractMultivariateProcess` and provide +implementations for the methods below. +- `mark_distribution(pp, t, h, d)` +- `ground_intensity(pp, t, h, d)` +- `integrated_ground_intensity(pp, t, h, d)` +- `ground_intensity_bound(pp, t, h, d)` +- `fit(Type{pp}, h)` +- `simulate(pp, h)` +In all the cases, `pp` is the point process being implemented, `t` is the instant in +which the function will be evaluated, `h` is the history up to `t` and `d` is the dimension. +Other methods should be implemented if permance is an issue. +The process is expected to have a field `mark_dist::Vector{<:AbstractMarkDistribution}`. If this +field is not present, `Base.ndims(pp)` must also be implemented. """ abstract type AbstractMultivariateProcess <: AbstractPointProcess end Base.ndims(pp::AbstractMultivariateProcess) = length(pp.mark_dist) -function mark_distribution(pp::AbstractMultivariateProcess, t, h) +function mark_distribution(pp::AbstractMultivariateProcess, t, h::History) return [mark_distribution(pp, t, h, d) for d in 1:ndims(pp)] end @@ -23,7 +37,7 @@ function ground_intensity_bound(pp::AbstractMultivariateProcess, t, h::History) return [ground_intensity_bound(pp, t, h, d) for d in 1:ndims(pp)] end -function intensity(pp::AbstractMultivariateProcess, m, t, h, d) +function intensity(pp::AbstractMultivariateProcess, m, t, h::History, d::Int) return ground_intensity(pp, t, h, d) * densityof(pp.mark_dist[d], t, h, m) end @@ -31,10 +45,26 @@ function intensity(pp::AbstractMultivariateProcess, m, t, h::History) return [intensity(pp, m, t, h, d) for d in 1:ndims(pp)] end -function log_intensity(pp::AbstractMultivariateProcess, m, t, h, d) +function log_intensity(pp::AbstractMultivariateProcess, m, t, h::History, d::Int) return log(intensity(pp, m, t, h, d)) end function log_intensity(pp::AbstractMultivariateProcess, m, t, h::History) return [log_intensity(pp, m, t, h, d) for d in 1:ndims(pp)] end + +function time_change(h::History{T}, pp::AbstractMultivariateProcess) where {T} + tmin = typemax(T) + tmax = typemin(T) + transformed_times = [zeros(T, nb_events(h, d)) for d in 1:ndims(h)] + for d in 1:ndims(pp) + Λ(t) = integrated_ground_intensity(pp, h, min_time(h), t, d) + transformed_times[d] .= Λ.(event_times(h, d)) + new_tmin = Λ(min_time(h)) + tmin = new_tmin < tmin ? new_tmin : tmin + new_tmax = Λ(max_time(h)) + tmax = new_tmax > tmax ? new_tmax : tmax + end + transformed_marks = [collect(event_marks(h, d)) for d in 1:ndims(h)] + return History(transformed_times, tmin, tmax, transformed_marks) +end diff --git a/src/simulation.jl b/src/simulation.jl index dfed8ca..e98d9f8 100644 --- a/src/simulation.jl +++ b/src/simulation.jl @@ -42,6 +42,10 @@ function simulate_ogata( return h end +function simulate(rng::AbstractRNG, pp::AbstractUnivariateProcess, args...; kwargs...) + return simulate_ogata(rng, pp, args...; kwargs...) +end + function simulate(pp::AbstractPointProcess, args...; kwargs...) return simulate(default_rng(), pp, args...; kwargs...) end diff --git a/src/univariate/hawkes/fit.jl b/src/univariate/hawkes/fit.jl index c4ec8b6..bd67280 100644 --- a/src/univariate/hawkes/fit.jl +++ b/src/univariate/hawkes/fit.jl @@ -1,5 +1,5 @@ """ - StatsAPI.fit(rng, ::Type{HawkesProcess{T}}, h::History; step_tol::Float64 = 1e-6, max_iter::Int = 1000) where {T<:Real} + StatsAPI.fit(::Type{HawkesProcess{T}}, h::History; step_tol::Float64 = 1e-6, max_iter::Int = 1000, rng::AbstractRNG=default_rng()) where {T<:Real} Expectation-Maximization algorithm from [Lewis2011](@cite). The relevant calculations are in page 4, equations 6-13. @@ -22,8 +22,7 @@ The algorithm consists in the following steps: ω = D / div 5. If convergence criterion is met, return updated parameters, otherwise, back to step 2. -Notice that, in the implementation, the process is normalized so the average inter-event time is equal to 1 and, -therefore, the interval of the process is transformed from T to N. Also, in equation (8) in the paper, +Notice that, in the implementation, the process is normalized so the average inter-event time is equal to 1 and, therefore, the interval of the process is transformed from T to N. Also, in equation (8) in the paper, ∑_{i=1:n} pᵢᵢ = ∑_{i=1:n} (1 - ∑_{j < i} Dᵢⱼ) = N - D. @@ -35,8 +34,10 @@ function StatsAPI.fit( max_iter::Int=1000, rng::AbstractRNG=default_rng(), ) where {T<:Real,MD<:PointProcessMarkDistribution} + mark_dist = fit(MD, h) + n = nb_events(h) - n == 0 && return HawkesProcess(zero(T), zero(T), zero(T)) + n == 0 && return HawkesProcess(zero(T), zero(T), zero(T), mark_dist) tmax = T(duration(h)) # Normalize times so average inter-event time is 1 (T -> n) @@ -95,7 +96,6 @@ function StatsAPI.fit( n_iters >= max_iter && @warn "Maximum number of iterations reached without convergence." - mark_dist = fit(MD, h) # Unnormalize back to original time scale (T -> tmax): # parameters in normalized space (') relate to original by μ0=μ'*(n/tmax), ω0=ω'*(n/tmax), α0=(ψ'ω')*(n/tmax) return HawkesProcess(μ * (n / tmax), ψ * ω * (n / tmax), ω * (n / tmax), mark_dist) diff --git a/src/univariate/hawkes/hawkes_process.jl b/src/univariate/hawkes/hawkes_process.jl index f806862..1d0269d 100644 --- a/src/univariate/hawkes/hawkes_process.jl +++ b/src/univariate/hawkes/hawkes_process.jl @@ -43,22 +43,22 @@ end HawkesProcess(μ, α, ω) = HawkesProcess(μ, α, ω, NoMarks()) -function ground_intensity(hp::HawkesProcess, h::History, t) +function ground_intensity(hp::HawkesProcess, t, h::History) activation = sum(exp.(hp.ω .* (@view h.times[1:(searchsortedfirst(h.times, t) - 1)]))) return hp.μ + (hp.α * activation / exp(hp.ω * t)) end -function integrated_ground_intensity(hp::HawkesProcess{T}, h::History, tmin, tmax) where {T} - U = promote_type(T, typeof(tmin), typeof(tmax)) - times = event_times(h, h.tmin, tmax) +function integrated_ground_intensity(hp::HawkesProcess{T}, h::History, a, b) where {T} + U = promote_type(T, typeof(a), typeof(b)) + times = event_times(h, h.tmin, b) integral = zero(U) for ti in times # Integral of activation function. 'max(tmin - ti, 0)' corrects for events that occurred # inside or outside the interval [tmin, tmax]. - integral += (exp(-hp.ω * max(tmin - ti, 0)) - exp(-hp.ω * (tmax - ti))) + integral += (exp(-hp.ω * max(a - ti, 0)) - exp(-hp.ω * (b - ti))) end integral *= hp.α / hp.ω - integral += hp.μ * (tmax - tmin) # Integral of base rate + integral += hp.μ * (b - a) # Integral of base rate return integral end @@ -70,7 +70,7 @@ function DensityInterface.logdensityof(hp::HawkesProcess, h::History) return sum(log.(hp.μ .+ (hp.α .* A))) - # Value of intensity at each event (hp.μ * duration(h)) - # Integral of base rate ((hp.α / hp.ω) * sum(1 .- exp.(-hp.ω .* (duration(h) .- h.times)))) + # Integral of each kernel - sum(log.([densityof(hp.mark_dist, t, h, m) for (t, m) in zip(h.times, h.marks)])) # Lok likelihood of marks + sum(log.([densityof(hp.mark_dist, t, h, m) for (t, m) in zip(h.times, h.marks)])) # Loglikelihood of marks end function time_change(h::History{R,M}, hp::HawkesProcess) where {R<:Real,M} diff --git a/src/univariate/hawkes/simulation.jl b/src/univariate/hawkes/simulation.jl index f436e92..f70afd8 100644 --- a/src/univariate/hawkes/simulation.jl +++ b/src/univariate/hawkes/simulation.jl @@ -3,9 +3,11 @@ function simulate(rng::AbstractRNG, hp::HawkesProcess, tmin, tmax) sim_desc = generate_descendants(rng, sim, tmax, hp.α, hp.ω) # Recursively generates descendants from first events append!(sim, sim_desc) sort!(sim) - h_temp = History(sim, tmin, tmax) - marks = [sample_mark(hp.mark_dist, t, h_temp) for t in event_times(h_temp)] - return History(; times=sim, tmin=tmin, tmax=tmax, marks=marks) + h = History(eltype(sim)[], tmin, tmax, eltype(hp.mark_dist)[], Nothing[], 1) + for event in sim + push!(h, event, sample_mark(hp, event, h)) + end + return h end #= diff --git a/src/univariate/poisson/fit.jl b/src/univariate/poisson/fit.jl index 52ade53..1825128 100644 --- a/src/univariate/poisson/fit.jl +++ b/src/univariate/poisson/fit.jl @@ -1,4 +1,9 @@ ## Fit MLE +function StatsAPI.fit(::Type{PoissonProcess{R,D}}, h::History; kwargs...) where {R,D} + mark_dist = fit(D, h.marks) + λ = nb_events(h) / duration(h) + return PoissonProcess(R(λ), mark_dist) +end function StatsAPI.fit( ::Type{PoissonProcess{R,D}}, ss::PoissonProcessStats{R1,R2}; kwargs... @@ -8,8 +13,10 @@ function StatsAPI.fit( return PoissonProcess(λ, mark_dist) end -function StatsAPI.fit(pptype::Type{PoissonProcess{R,D}}, args...; kwargs...) where {R,D} - ss = suffstats(pptype, args...) +function StatsAPI.fit( + pptype::Type{PoissonProcess{R,D}}, hs::AbstractVector{<:History}; kwargs... +) where {R,D} + ss = suffstats(pptype, hs) return fit(pptype, ss) end diff --git a/src/univariate/poisson/inhomogeneous/inhomogeneous_poisson_process.jl b/src/univariate/poisson/inhomogeneous/inhomogeneous_poisson_process.jl index e922db4..58607b1 100644 --- a/src/univariate/poisson/inhomogeneous/inhomogeneous_poisson_process.jl +++ b/src/univariate/poisson/inhomogeneous/inhomogeneous_poisson_process.jl @@ -61,9 +61,9 @@ end ## AbstractPointProcess interface -ground_intensity(pp::InhomogeneousPoissonProcess, t, h) = pp.intensity_function(t) +ground_intensity(pp::InhomogeneousPoissonProcess, t, h::History) = pp.intensity_function(t) -function intensity(pp::InhomogeneousPoissonProcess, m, t, h) +function intensity(pp::InhomogeneousPoissonProcess, m, t, h::History) return ground_intensity(pp, t, h) * densityof(pp.mark_dist, t, h, m) end diff --git a/src/univariate_process.jl b/src/univariate_process.jl index fd91d85..4e8c0fd 100644 --- a/src/univariate_process.jl +++ b/src/univariate_process.jl @@ -7,18 +7,32 @@ abstract type AbstractUnivariateProcess <: AbstractPointProcess end Base.ndims(::AbstractUnivariateProcess) = 1 -function mark_distribution(pp::AbstractUnivariateProcess, t, h) +function mark_distribution(pp::AbstractUnivariateProcess, t, h::History) return mark_distribution(pp.mark_dist, t, h) end -function intensity(pp::AbstractUnivariateProcess, m, t, h) +function sample_mark(rng::AbstractRNG, pp::AbstractPointProcess, t, h::History) + return sample_mark(rng, pp.mark_dist, t, h) +end +sample_mark(pp::AbstractPointProcess, t, h::History) = sample_mark(default_rng(), pp, t, h) + +function intensity(pp::AbstractUnivariateProcess, m, t, h::History) return ground_intensity(pp, t, h) * densityof(pp.mark_dist, t, h, m) end -function log_intensity(pp::AbstractUnivariateProcess, m, t, h) +function log_intensity(pp::AbstractUnivariateProcess, m, t, h::History) return log(intensity(pp, m, t, h)) end +""" + logdensityof(pp, h) + +Compute the log probability density function for a temporal point process `pp` applied to history `h`: +``` +ℓ(h) = Σₖ log λ(tₖ|hₖ) - Λ(h) +``` +The default method uses a loop over events combined with `integrated_ground_intensity`, but it should be reimplemented for specific processes if faster computation is possible. +""" function DensityInterface.logdensityof(pp::AbstractUnivariateProcess, h::History) l = -integrated_ground_intensity(pp, h, min_time(h), max_time(h)) for (t, m) in zip(event_times(h), event_marks(h)) @@ -27,7 +41,7 @@ function DensityInterface.logdensityof(pp::AbstractUnivariateProcess, h::History return l end -function time_change(h::History, pp::AbstractUnivariateProcess) +function time_change(h::History{T}, pp::AbstractUnivariateProcess) where {T} Λ(t) = integrated_ground_intensity(pp, h, min_time(h), t) return time_change(h, Λ) end diff --git a/test/bounded_point_process.jl b/test/bounded_point_process.jl index 718b9ae..e5ece17 100644 --- a/test/bounded_point_process.jl +++ b/test/bounded_point_process.jl @@ -14,20 +14,3 @@ h = simulate(rng, bpp) @test all(ground_intensity_bound(bpp, 243.0, h, 1) .≈ (intensities[1], Inf)) @test integrated_ground_intensity(bpp, h, 342, 598) ≈ intensities * (598 - 342) - -struct FakePoisson <: AbstractUnivariateProcess - λ::Float64 - mark_dist::NoMarks -end - -PointProcesses.ground_intensity(fp::FakePoisson, t, h::History) = fp.λ -function PointProcesses.integrated_ground_intensity(fp::FakePoisson, h::History, a, b) - return fp.λ * (b - a) -end -PointProcesses.ground_intensity_bound(fp::FakePoisson, t, h::History) = (fp.λ, Inf) - -pp = BoundedPointProcess(FakePoisson(1.0, NoMarks()), 0.0, 10.0) -h = History([1.0, 2.0], 0.0, 10.0) - -@test simulate(pp) isa History -@test time_change(h, pp).times == h.times diff --git a/test/custom_processes.jl b/test/custom_processes.jl new file mode 100644 index 0000000..d6813d8 --- /dev/null +++ b/test/custom_processes.jl @@ -0,0 +1,72 @@ +@testset "Univariate" begin + struct FakePoisson <: AbstractUnivariateProcess + λ::Float64 + mark_dist::NoMarks + end + + PointProcesses.ground_intensity(fp::FakePoisson, t, h::History) = fp.λ + + function PointProcesses.integrated_ground_intensity(fp::FakePoisson, h::History, a, b) + return fp.λ * (b - a) + end + + PointProcesses.ground_intensity_bound(fp::FakePoisson, t, h::History) = (fp.λ, Inf) + + fp = FakePoisson(1.0, NoMarks()) + h = simulate(fp, 0.0, 10.0) + + @test h isa History + @test ndims(h) == 1 + @test ndims(fp) == 1 + @test intensity(fp, nothing, 1.0, h) == 1.0 + @test log_intensity(fp, nothing, 1.0, h) == 0.0 + @test mark_distribution(fp, 1.0, h) == Dirac(nothing) + @test sample_mark(fp, 1.0, h) === nothing + @test time_change(h, fp).times == h.times + @test logdensityof(fp, h) == -10.0 + + fp_bounded = BoundedPointProcess(FakePoisson(1.0, NoMarks()), 0.0, 10.0) + @test simulate(fp_bounded) isa History + @test time_change(h, fp_bounded).times == h.times +end + +@testset "Multivariate" begin + struct FakeMultivariatePoisson <: AbstractMultivariateProcess + λ::Vector{Float64} + mark_dist::Vector{NoMarks} + end + + PointProcesses.mark_distribution(::FakeMultivariatePoisson, t, ::History, ::Int) = + Dirac(nothing) + + PointProcesses.ground_intensity(fmp::FakeMultivariatePoisson, t, h::History, d::Int) = + fmp.λ[d] + + function PointProcesses.integrated_ground_intensity( + fmp::FakeMultivariatePoisson, h::History, a, b, d::Int + ) + return fmp.λ[d] * (b - a) + end + + PointProcesses.ground_intensity_bound( + fmp::FakeMultivariatePoisson, t, h::History, d::Int + ) = (fmp.λ[d], Inf) + + fmp = FakeMultivariatePoisson([1.0, 2.0], [NoMarks(), NoMarks()]) + h = History([rand(10), rand(10)], 0.0, 1.0) + + @test h isa History + @test ndims(h) == 2 + @test ndims(fmp) == 2 + @test mark_distribution(fmp, 1.0, h) == [Dirac(nothing), Dirac(nothing)] + @test ground_intensity(fmp, 1.0, h) == [1.0, 2.0] + @test integrated_ground_intensity(fmp, h, 0.0, 2.0) == [2.0, 4.0] + @test ground_intensity_bound(fmp, 1.0, h) == [(1.0, Inf), (2.0, Inf)] + @test intensity(fmp, nothing, 1.0, h, 1) == 1.0 + @test intensity(fmp, nothing, 1.0, h) == [1.0, 2.0] + @test log_intensity(fmp, nothing, 1.0, h, 1) == 0.0 + @test log_intensity(fmp, nothing, 1.0, h) == [0.0, log(2.0)] + h_transf = time_change(h, fmp) + @test event_times(h_transf, 1) == event_times(h, 1) + @test event_times(h_transf, 2) == event_times(h, 2) .* 2 +end diff --git a/test/hawkes.jl b/test/hawkes.jl index 65a277a..1b3a609 100644 --- a/test/hawkes.jl +++ b/test/hawkes.jl @@ -23,8 +23,8 @@ integral = @test isa(time_change(h_big, hp), typeof(h_big)) # Ground intensity -@test ground_intensity(hp, h, 1) == 1 -@test ground_intensity(hp, h, 2) == 1 + hp.α * exp(-hp.ω * 1) +@test ground_intensity(hp, 1, h) == 1 +@test ground_intensity(hp, 2, h) == 1 + hp.α * exp(-hp.ω * 1) # Integrated ground intensity @test integrated_ground_intensity(hp, h, h.tmin, h.tmax) ≈ integral diff --git a/test/runtests.jl b/test/runtests.jl index 4c0d52e..1ac0906 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -50,6 +50,9 @@ DocMeta.setdocmeta!(PointProcesses, :DocTestSetup, :(using PointProcesses); recu @testset verbose = true "Bounded" begin include("bounded_point_process.jl") end + @testset verbose = true "CustomProcesses" begin + include("custom_processes.jl") + end @testset verbose = true "IndependentMultivariate" begin include("independent_multivariate.jl") end