From e98f3e6cc5b4111866c973f05d1a2d88238a0999 Mon Sep 17 00:00:00 2001 From: jkling Date: Thu, 28 May 2026 19:44:11 +0200 Subject: [PATCH 1/6] 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/6] 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/6] 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/6] 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/6] 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 From 5d08bc99540aacb47596ee2e15f8d1daf7c88d10 Mon Sep 17 00:00:00 2001 From: jkling Date: Sat, 18 Jul 2026 11:59:22 +0200 Subject: [PATCH 6/6] Fix for PR - Fixed constructor for empty history - History does not deduce number of dimensions - Small fixes for Hawkes processes --- docs/examples/CustomProcesses.jl | 3 +- src/bounded_point_process.jl | 7 ++ src/history.jl | 99 ++++++++++++------- src/mark_distributions.jl | 2 +- src/multivariate_process.jl | 10 ++ src/simulation.jl | 2 +- src/univariate/hawkes/hawkes_process.jl | 23 +++-- src/univariate/hawkes/simulation.jl | 2 +- .../inhomogeneous_poisson_process.jl | 4 - src/univariate/poisson/simulation.jl | 4 +- src/univariate_process.jl | 3 +- test/history.jl | 35 +++++-- test/inhomogeneous_poisson_process.jl | 14 +-- test/mark_distributions.jl | 2 +- 14 files changed, 139 insertions(+), 71 deletions(-) diff --git a/docs/examples/CustomProcesses.jl b/docs/examples/CustomProcesses.jl index a7a96c1..c3e9312 100644 --- a/docs/examples/CustomProcesses.jl +++ b/docs/examples/CustomProcesses.jl @@ -173,7 +173,8 @@ function StatsAPI.fit( optimal = Optim.minimizer(result) solved = true - catch DomainError + catch e + e isa DomainError || rethrow() n_tries += 1 end end diff --git a/src/bounded_point_process.jl b/src/bounded_point_process.jl index c123318..c2b236e 100644 --- a/src/bounded_point_process.jl +++ b/src/bounded_point_process.jl @@ -40,6 +40,13 @@ function integrated_ground_intensity(bpp::BoundedPointProcess, args...) return integrated_ground_intensity(bpp.pp, args...) end +function DensityInterface.logdensityof(bpp::BoundedPointProcess, h) + if min_time(h) < min_time(bpp) || max_time(h) > max_time(bpp) + throw(ArgumentError("History is defined outside the bounds of the process")) + end + return logdensityof(bpp.pp, h) +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")) diff --git a/src/history.jl b/src/history.jl index d0ef40e..297ab9f 100644 --- a/src/history.jl +++ b/src/history.jl @@ -41,9 +41,9 @@ struct History{T<:Real,M,D} times::AbstractVector{R1}, tmin::R2, tmax::R3, - marks::AbstractVector{M}=fill(nothing, length(times)), - dims::AbstractVector{D}=fill(nothing, length(times)), - N::Int=length(unique(dims)); + marks::AbstractVector{M}, + dims::AbstractVector{D}, + N::Int; check_args=true, ) where {R1<:Real,R2<:Real,R3<:Real,M,D} if check_args @@ -146,42 +146,72 @@ end function History( times::AbstractVector{<:AbstractVector{R1}}, tmin::R2, tmax::R3; check_args=true ) where {R1<:Real,R2<:Real,R3<:Real} - nots = [fill(nothing, length(times[i])) for i in 1:length(times)] - return History(times, tmin, tmax, nots; check_args=check_args) + marks = [fill(nothing, length(times[i])) for i in eachindex(times)] + return History(times, tmin, tmax, marks; check_args=check_args) 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))) - if dims === nothing - dims = fill(nothing, length(times)) - N = 1 - else - N = length(unique(dims)) +function History( + times::AbstractVector{R1}, tmin::R2, tmax::R3, marks::AbstractVector; check_args=true +) where {R1<:Real,R2<:Real,R3<:Real} + dims = fill(nothing, length(times)) + return History(times, tmin, tmax, marks, dims, 1; check_args=check_args) +end + +function History( + times::AbstractVector{R1}, tmin::R2, tmax::R3; check_args=true +) where {R1<:Real,R2<:Real,R3<:Real} + marks = fill(nothing, length(times)) + return History(times, tmin, tmax, marks; check_args=check_args) +end + +function History(; times, tmin, tmax, marks=nothing, check_args=true) + if times isa AbstractVector{<:Real} + if isnothing(marks) + marks = fill(nothing, length(times)) end - return History(times, tmin, tmax, marks, dims, N; check_args=check_args) + dims = fill(nothing, length(times)) + return History(times, tmin, tmax, marks, dims, 1; check_args=check_args) else - marks === nothing && - (marks = [fill(nothing, length(times[i])) for i in 1:length(times)]) + if isnothing(marks) + marks = [fill(nothing, length(times[i])) for i in eachindex(times)] + end return History(times, tmin, tmax, marks; check_args=check_args) end end -function History(tmin::R1, tmax::R2, N::Int=1) where {R1<:Real,R2<:Real} +function History(tmin::R1, tmax::R2, M::Type) where {R1,R2<:Real} R = promote_type(R1, R2) - return History(R[], tmin, tmax, [], [], N) + return History(R[], tmin, tmax, M[], Nothing[], 1; check_args=false) +end + +function History(tmin::R1, tmax::R2, M::Type, N::Int) where {R1,R2<:Real} + if N == 1 + return History(tmin, tmax, M) + else + R = promote_type(R1, R2) + return History(R[], tmin, tmax, M[], Int[], N; check_args=false) + end end function History(h::History, d::Int) times = event_times(h, d) marks = event_marks(h, d) - return History(times, min_time(h), max_time(h), marks) + dims = fill(nothing, length(times)) + return History(times, min_time(h), max_time(h), marks, dims, 1; check_args=false) end function Base.show(io::IO, h::History{T,M}) where {T,M} - return print( - io, "History{$T,$M} with $(nb_events(h)) events on interval [$(h.tmin), $(h.tmax))" - ) + if ndims(h) == 1 + return print( + io, + "Univariate History{$T,$M} with $(nb_events(h)) events on interval [$(h.tmin), $(h.tmax))", + ) + else + return print( + io, + "$(ndims(h))-dimensional History{$T,$M} with $(nb_events(h)) events on interval [$(h.tmin), $(h.tmax))", + ) + end end """ @@ -455,18 +485,21 @@ h1 coincides with the beginning of h2, then create a new event history by concatenating h1 and h2. """ function Base.cat(h1::History, h2::History) - max_time(h1) ≈ min_time(h2) || throw( - DomainError( - (h1.tmax, h2.tmin), - "End of h1's interval must coincide with start of h2's interval", - ), - ) + if !(max_time(h1) ≈ min_time(h2)) + throw( + DomainError( + (h1.tmax, h2.tmin), + "End of h1's interval must coincide with start of h2's interval", + ), + ) + end + if ndims(h1) != ndims(h2) + throw(DimensionMismatch("h1 and h2 must have the same number of dimensions")) + end times = [h1.times; h2.times] marks = [h1.marks; h2.marks] dims = [h1.dims; h2.dims] - return History(; - times=times, tmin=h1.tmin, tmax=h2.tmax, marks=marks, dims=dims, check_args=false - ) + return History(times, h1.tmin, h2.tmax, marks, dims, ndims(h1); check_args=false) end function time_change(h::History{T}, Λ) where {T} @@ -492,9 +525,7 @@ function split_into_chunks(h::History{T,M,D}, chunk_duration) where {T,M,D} times = [t for t in event_times(h) if a <= t < b] marks = [m for (t, m) in zip(event_times(h), event_marks(h)) if a <= t < b] dims = [d for (t, d) in zip(event_times(h), event_dims(h)) if a <= t < b] - chunk = History(; - times=times, marks=marks, dims=dims, tmin=a, tmax=b, check_args=false - ) + chunk = History(times, a, b, marks, dims, h.N; check_args=false) push!(chunks, chunk) end return chunks diff --git a/src/mark_distributions.jl b/src/mark_distributions.jl index 8174fcd..74aa06d 100644 --- a/src/mark_distributions.jl +++ b/src/mark_distributions.jl @@ -36,7 +36,7 @@ end "The type of the marks returned by the mark distribution" function Base.eltype(md::AbstractMarkDistribution) - return typeof(sample_mark(md, 0.0, History(0.0, 1.0))) + return typeof(sample_mark(md, 0.0, History(0.0, 1.0, Nothing))) end "The likelihood of a mark `m` occurring in an event at time `t` after history `h`" diff --git a/src/multivariate_process.jl b/src/multivariate_process.jl index a6425b7..9d9b6c6 100644 --- a/src/multivariate_process.jl +++ b/src/multivariate_process.jl @@ -25,6 +25,16 @@ function mark_distribution(pp::AbstractMultivariateProcess, t, h::History) return [mark_distribution(pp, t, h, d) for d in 1:ndims(pp)] end +function sample_mark( + rng::AbstractRNG, pp::AbstractMultivariateProcess, t, h::History, d::Int +) + return sample_mark(rng, pp.mark_dist[d], t, h) +end + +function sample_mark(pp::AbstractMultivariateProcess, t, h::History, d::Int) + return sample_mark(default_rng(), pp, t, h, d) +end + function ground_intensity(pp::AbstractMultivariateProcess, t, h::History) return [ground_intensity(pp, t, h, d) for d in 1:ndims(pp)] end diff --git a/src/simulation.jl b/src/simulation.jl index e98d9f8..144b308 100644 --- a/src/simulation.jl +++ b/src/simulation.jl @@ -20,7 +20,7 @@ function simulate_ogata( rng::AbstractRNG, pp::AbstractUnivariateProcess, tmin::T, tmax::T ) where {T<:Real} M = eltype(pp.mark_dist) - h = History(; times=T[], marks=M[], tmin=tmin, tmax=tmax) + h = History(tmin, tmax, M) t = tmin while t < tmax B, L = ground_intensity_bound(pp, t + eps(t), h) diff --git a/src/univariate/hawkes/hawkes_process.jl b/src/univariate/hawkes/hawkes_process.jl index 1d0269d..5a2be11 100644 --- a/src/univariate/hawkes/hawkes_process.jl +++ b/src/univariate/hawkes/hawkes_process.jl @@ -31,10 +31,11 @@ struct HawkesProcess{T<:Real,D<:PointProcessMarkDistribution} <: AbstractUnivari function HawkesProcess( μ::T1, α::T2, ω::T3, mark_dist::D ) where {T1,T2,T3<:Real,D<:PointProcessMarkDistribution} - any((μ, α, ω) .< 0) && + if any((μ, α, ω) .< 0) throw(DomainError((μ, α, ω), "All parameters must be non-negative.")) - (α > 0 && α >= ω) && - throw(DomainError((α, ω), "Parameter ω must be strictly smaller than α")) + elseif (α > 0 && α >= ω) + throw(DomainError((α, ω), "Parameter α must be strictly smaller than ω")) + end T = promote_type(T1, T2, T3) (μ_T, α_T, ω_T) = convert.(T, (μ, α, ω)) return new{T,D}(μ_T, α_T, ω_T, mark_dist) @@ -44,8 +45,9 @@ end HawkesProcess(μ, α, ω) = HawkesProcess(μ, α, ω, NoMarks()) 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)) + past_times = event_times(h, h.tmin, t) + activation = sum(exp.(-hp.ω .* (t .- past_times))) + return hp.μ + hp.α * activation end function integrated_ground_intensity(hp::HawkesProcess{T}, h::History, a, b) where {T} @@ -67,10 +69,13 @@ function DensityInterface.logdensityof(hp::HawkesProcess, h::History) 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)])) # Loglikelihood of marks + sum_intensities_t = sum(log.(hp.μ .+ (hp.α .* A))) # Value of intensity at each event + integral_base = hp.μ * duration(h) + integral_kernels = (hp.α / hp.ω) * sum(1 .- exp.(-hp.ω .* (h.tmax .- h.times))) + loglike_marks = sum( + log(densityof(hp.mark_dist, t, h, m)) for (t, m) in zip(h.times, h.marks); init=0.0 + ) + return sum_intensities_t - integral_base - integral_kernels + loglike_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 f70afd8..bcde023 100644 --- a/src/univariate/hawkes/simulation.jl +++ b/src/univariate/hawkes/simulation.jl @@ -5,7 +5,7 @@ function simulate(rng::AbstractRNG, hp::HawkesProcess, tmin, tmax) sort!(sim) h = History(eltype(sim)[], tmin, tmax, eltype(hp.mark_dist)[], Nothing[], 1) for event in sim - push!(h, event, sample_mark(hp, event, h)) + push!(h, event, sample_mark(rng, hp, event, h)) end return h end diff --git a/src/univariate/poisson/inhomogeneous/inhomogeneous_poisson_process.jl b/src/univariate/poisson/inhomogeneous/inhomogeneous_poisson_process.jl index 58607b1..84884d5 100644 --- a/src/univariate/poisson/inhomogeneous/inhomogeneous_poisson_process.jl +++ b/src/univariate/poisson/inhomogeneous/inhomogeneous_poisson_process.jl @@ -63,10 +63,6 @@ end ground_intensity(pp::InhomogeneousPoissonProcess, t, h::History) = pp.intensity_function(t) -function intensity(pp::InhomogeneousPoissonProcess, m, t, h::History) - return ground_intensity(pp, t, h) * densityof(pp.mark_dist, t, h, m) -end - """ ground_intensity_bound(pp::InhomogeneousPoissonProcess, t, h) diff --git a/src/univariate/poisson/simulation.jl b/src/univariate/poisson/simulation.jl index bf4872e..a7a33ce 100644 --- a/src/univariate/poisson/simulation.jl +++ b/src/univariate/poisson/simulation.jl @@ -1,9 +1,9 @@ function simulate(rng::AbstractRNG, pp::PoissonProcess, tmin::T, tmax::T) where {T<:Real} - h = History(T[], tmin, tmax, eltype(pp.mark_dist)[]) + h = History(tmin, tmax, eltype(pp.mark_dist)) inter_dist = Exponential(inv(pp.λ)) t = T(rand(rng, inter_dist)) while t < tmax - m = sample_mark(pp.mark_dist, t, h) + m = sample_mark(pp, t, h) push!(h, t, m; check_args=false) t += rand(rng, inter_dist) end diff --git a/src/univariate_process.jl b/src/univariate_process.jl index 4e8c0fd..b6d062d 100644 --- a/src/univariate_process.jl +++ b/src/univariate_process.jl @@ -11,9 +11,10 @@ function mark_distribution(pp::AbstractUnivariateProcess, t, h::History) return mark_distribution(pp.mark_dist, t, h) end -function sample_mark(rng::AbstractRNG, pp::AbstractPointProcess, t, h::History) +function sample_mark(rng::AbstractRNG, pp::AbstractUnivariateProcess, 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) diff --git a/test/history.jl b/test/history.jl index 339dc5e..b4e49fb 100644 --- a/test/history.jl +++ b/test/history.jl @@ -1,15 +1,18 @@ @testset "Univariate History" begin # Constructors - h_empty1 = History(0.0, 1.0) - h_empty2 = History(0.0, 1.0, 2) + h_empty1 = History(0.0, 1.0, Nothing) + h_empty2 = History(0.0, 1.0, Nothing, 2) - @test h_empty1 isa History{Float64,Any} - @test h_empty2 isa History{Float64,Any} + @test h_empty1 isa History{Float64,Nothing} + @test h_empty2 isa History{Float64,Nothing} @test ndims(h_empty1) == 1 @test ndims(h_empty2) == 2 @test isempty(h_empty1) @test isempty(h_empty2) + @test History(; times=sort(rand(3)), tmin=0.0, tmax=1.0) isa History{Float64,Nothing} + @test History(; times=[sort(rand(3)), sort(rand(5))], tmin=0.0, tmax=1.0) isa + History{Float64,Nothing} h = History([0.2, 0.8, 1.1], 0.0, 2.0, ["a", "b", "c"]) @test h isa History{Float64,String} @@ -42,13 +45,15 @@ h2 = History(; times=[2.3], marks=["e"], tmin=2.0, tmax=2.5) - @test string(h2) == "History{Float64,String} with 1 events on interval [2.0, 2.5)" + @test string(h2) == + "Univariate History{Float64,String} with 1 events on interval [2.0, 2.5)" h_cat = cat(h, h2) @test nb_events(h_cat) == 5 @test duration(h_cat) == 2.5 @test length(split_into_chunks(h_cat, 0.3)) == 9 + @test_throws DomainError cat(h, h) h_exp = time_change(h_cat, exp) @@ -59,9 +64,9 @@ @test h2.times == [2.3, 2.4, 2.45] @test h2.marks == ["e", "f", "g"] - @test isa(History(rand(3), 0, BigFloat(1)), History{BigFloat,Nothing}) - @test_throws DomainError History(rand(3), 1, 0) - @test_throws DimensionMismatch History(rand(3), 0, 1, ["a", "b"]) + @test isa(History(sort(rand(3)), 0, BigFloat(1)), History{BigFloat,Nothing}) + @test_throws DomainError History(sort(rand(3)), 1, 0) + @test_throws DimensionMismatch History(sort(rand(3)), 0, 1, ["a", "b"]) # Strict validation: events outside [tmin, tmax) must throw, not silently # discard. Same for unsorted times and same-time-same-dim repeats. @test_throws DomainError History([0.1, 1.1], 0, 1) @@ -90,7 +95,7 @@ end @test event_marks(h_multi1) == marks1 @test event_marks(h_multi2) == marks2 - @test_throws DomainError History(rand(3), 0, 1, rand(3), [1, 2, 3], 2) + @test_throws DomainError History(sort(rand(3)), 0, 1, rand(3), [1, 2, 3], 2) @test event_dims(History([[0.5]], 0, 1)) == [nothing] @test_throws DomainError History( @@ -114,6 +119,18 @@ end @test nb_events(h_multi, 0.0, 0.3, 1) == 1 @test event_dims(h_multi, 0.0, 0.3) == [1, 2] + # Test cat with multiple dimension + h_multi_slide = History([times1 .+ 1.0, times2 .+ 1.0], 1.0, 2.0, [marks1, marks2]) + h_multi_cat = cat(h_multi, h_multi_slide) + @test ndims(h_multi_cat) == 2 + @test nb_events(h_multi_cat) == 8 + @test min_time(h_multi_cat) == 0.0 + @test max_time(h_multi_cat) == 2.0 + + h_uni = History(sort(rand(3)) .+ 1.0, 1.0, 2.0, ["a", "b", "c"]) + + @test_throws DimensionMismatch cat(h_multi, h_uni) + # Test push! with dimension @test_throws DomainError push!(h_multi, 1.0, "d", nothing) push!(h_multi, 0.85, "e", 1) diff --git a/test/inhomogeneous_poisson_process.jl b/test/inhomogeneous_poisson_process.jl index 0f9fc80..e673dee 100644 --- a/test/inhomogeneous_poisson_process.jl +++ b/test/inhomogeneous_poisson_process.jl @@ -973,7 +973,7 @@ end cov1 = t -> t cov2 = t -> sin(t) intensity_cov = LinearCovariateIntensity(1.0, [0.5, 0.2], [cov1, cov2]) - h_empty = History(Float64[], 0.0, 5.0, Float64[]) + h_empty = History(0.0, 5.0, Float64) pp = InhomogeneousPoissonProcess(intensity_cov, Normal()) # This should use numerical integration @@ -993,7 +993,7 @@ end @testset "Generic fallback for custom functions" begin # Test custom lambda function custom_func = t -> 2.0 + 0.5 * t^2 - h_empty = History(Float64[], 0.0, 4.0, Float64[]) + h_empty = History(0.0, 4.0, Float64) pp = InhomogeneousPoissonProcess(custom_func, Uniform()) # ∫₀⁴ (2 + 0.5*t²) dt = [2t + 0.5*t³/3]₀⁴ = 8 + 0.5*64/3 = 8 + 32/3 ≈ 18.667 @@ -1014,7 +1014,7 @@ end @testset "Intensity bounds edge cases" begin @testset "ExponentialIntensity bound with b > 0 (increasing)" begin intensity_inc = ExponentialIntensity(2.0, 0.1) - h_empty = History(Float64[], 0.0, 10.0, Float64[]) + h_empty = History(0.0, 10.0, Float64) pp = InhomogeneousPoissonProcess(intensity_inc, Normal()) B, L = ground_intensity_bound(pp, 0.0, h_empty) @@ -1029,7 +1029,7 @@ end @testset "ExponentialIntensity bound with b < 0 (decreasing)" begin intensity_dec = ExponentialIntensity(5.0, -0.2) - h_empty = History(Float64[], 0.0, 10.0, Float64[]) + h_empty = History(0.0, 10.0, Float64) pp = InhomogeneousPoissonProcess(intensity_dec, Normal()) B, L = ground_intensity_bound(pp, 2.0, h_empty) @@ -1044,7 +1044,7 @@ end @testset "ExponentialIntensity bound with b ≈ 0 (constant)" begin intensity_const = ExponentialIntensity(3.0, 1e-12) - h_empty = History(Float64[], 0.0, 10.0, Float64[]) + h_empty = History(0.0, 10.0, Float64) pp = InhomogeneousPoissonProcess(intensity_const, Normal()) B, L = ground_intensity_bound(pp, 5.0, h_empty) @@ -1062,7 +1062,7 @@ end intensity_func = PolynomialIntensity([2.0, 0.3]) mark_dist = Normal(0.0, 1.0) pp = InhomogeneousPoissonProcess(intensity_func, mark_dist) - h_empty = History(Float64[], 0.0, 10.0, Float64[]) + h_empty = History(0.0, 10.0, Float64) t = 5.0 m = 1.5 @@ -1086,7 +1086,7 @@ end intensity_func = ExponentialIntensity(3.0, 0.05) mark_dist = Categorical([0.3, 0.5, 0.2]) pp = InhomogeneousPoissonProcess(intensity_func, mark_dist) - h_empty = History(Float64[], 0.0, 10.0, Float64[]) + h_empty = History(0.0, 10.0, Float64) t = 2.0 # ground_intensity at t=2: λ(2) = 3*exp(0.05*2) = 3*exp(0.1) diff --git a/test/mark_distributions.jl b/test/mark_distributions.jl index f0abd21..897164d 100644 --- a/test/mark_distributions.jl +++ b/test/mark_distributions.jl @@ -30,7 +30,7 @@ end t::Float64 end - h = History(0.0, 1.0) + h = History(0.0, 1.0, Float64) md = TestMark(1.0) @test_throws "not implemented" mark_distribution(md, 0.0, h)