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,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)
-    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
 
@@ -168,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

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,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)

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,12 +469,7 @@ 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, Λ)
+function time_change(h::History{T}, Λ) where {T}
     new_times = Λ.(event_times(h))
     new_marks = copy(event_marks(h))
     new_tmin = Λ(min_time(h))

src/multivariate/independent_multivariate.jl

@@ -18,56 +18,28 @@ 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::Int)
     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)]
-end
-
-function mark_distribution(pp::IndependentMultivariateProcess, t, h, 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)
+function intensity(pp::IndependentMultivariateProcess, m, t, h::History, d::Int)
     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)
-    return log(intensity(pp, m, t, h, d))
+function mark_distribution(pp::IndependentMultivariateProcess, t, h::History, d::Int)
+    return mark_distribution(pp.processes[d], t, History(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::Int)
     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)
+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
 
-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(

src/multivariate_process.jl

@@ -0,0 +1,70 @@
+"""
+    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::History)
+    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::History, d::Int)
+    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::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