Interface for Multivariate processes

docs/examples/CustomProcesses.jl

@@ -33,8 +33,8 @@ using PointProcesses
 using Distributions
 using Plots
 using StatsAPI
+using StatsBase
 using Optim
-
 import PointProcesses: NoMarks, AbstractMarkDistribution, AbstractUnivariateProcess
 
 struct TwoStateModel{R<:Real,D<:PointProcessMarkDistribution} <: AbstractUnivariateProcess
@@ -141,15 +141,46 @@ 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)
-    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 e
+            e isa DomainError || rethrow()
+            n_tries += 1
+        end
+    end
+    if !solved
+        throw(ErrorException("Solver did not converge."))
+    end
     return TwoStateModel(optimal..., NoMarks())
 end
 
@@ -168,8 +199,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

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

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

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)

src/PointProcesses.jl

@@ -86,6 +86,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")
 
@@ -103,6 +105,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
 

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
 
@@ -113,23 +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.
-"""
-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
-
 """
     fit(::Type{PP}, h)
     fit(::Type{PP}, histories)
@@ -150,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)
 

src/bounded_point_process.jl

@@ -39,3 +39,19 @@ end
 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"))
+    end
+    return time_change(h, pp.pp)
+end
+
+simulate(pp::BoundedPointProcess) = simulate(pp.pp, pp.tmin, pp.tmax)

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,37 +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)))
-        dims === nothing && (dims = fill(nothing, length(times)))
-        return History(times, tmin, tmax, marks, dims; check_args=check_args)
+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
+        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
 
 """
@@ -450,26 +485,24 @@ 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
 
-"""
-    time_change(h, Λ)
-
-Apply the time rescaling `t -> Λ(t)` to history `h`.
-"""
-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))
@@ -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