Inhomogeneous point process

[rsenne]

Start of a PR for adding IPP models. Should resolve #33. I think this is a reasonably good start but before continuing, I think I'd like an extra set of eyes to make sure everyone is happy with the current structure. We can almsot certainly do better on the numerical integration so am open to suggestions. I also feel we may want to expand the AutoDiff support (maybe via DifferentationInterface).

This pull request introduces support for inhomogeneous Poisson processes (IPPs) to the point process modeling library. It adds new intensity function types, a flexible process type for IPPs, and comprehensive fitting methods using maximum likelihood estimation (MLE) with automatic differentiation. These changes enable modeling, simulation, and fitting of temporal Poisson processes with non-constant, parametric, or piecewise intensity functions.

Inhomogeneous Poisson Process Support

• Added the new InhomogeneousPoissonProcess type, representing Poisson processes with time-varying intensity and mark distributions, along with interface methods and documentation.
• Exported InhomogeneousPoissonProcess and several new intensity function types (PolynomialIntensity, ExponentialIntensity, SinusoidalIntensity, PiecewiseConstantIntensity, LinearCovariateIntensity) to the main package API.

Intensity Function Infrastructure

• Implemented parametric intensity function types for IPPs, including polynomial, exponential, sinusoidal, piecewise constant, and linear covariate forms, each with constructors, validation, and callable behavior.

Fitting and Estimation

• Added MLE-based fitting methods for IPPs, supporting automatic differentiation and numerical optimization via Optim.jl. Includes specialized fitting routines for each parametric intensity type and support for fitting from multiple event histories.

Integration and Simulation

• Integrated new IPP modules and intensity functions into the package, making them available for use and simulation.
• Added dependency on the Optim.jl package for numerical optimization in fitting routines.

[codecov]

Codecov Report

❌ Patch coverage is 98.32636% with 4 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
src/poisson/inhomogeneous/parametric_intensity.jl	0.00%	2 Missing ⚠️
src/poisson/inhomogeneous/fit.jl	97.67%	1 Missing ⚠️
src/poisson/inhomogeneous/intensity_methods.jl	98.96%	1 Missing ⚠️

📢 Thoughts on this report? Let us know!

[JoseKling]

I like the structure. Just to see if I get it.

If I want to define a point process with an arbitrary intensity function f(t), I just need to define, for example,

ipp = InhomogeneousPointProcess(f, Uniform())

and if I want to use the fit method, I just need to define the methods param_to_intensity and intensity_to_param and call the fit method with some initial guess. Is that it?

[rsenne]

I like the structure. Just to see if I get it.

If I want to define a point process with an arbitrary intensity function f(t), I just need to define, for example,

ipp = InhomogeneousPointProcess(f, Uniform())

and if I want to use the fit method, I just need to define the methods param_to_intensity and intensity_to_param and call the fit method with some initial guess. Is that it?

That's basically it with a small caveat--as you pointed out, assuming you have an arbitrary f(t), you can indeed call

ipp = InhomogeneousPointProcess(f, Uniform())

and you need to provide param_to_intensity and intensity_to_param. However, one would not call fitbut fit_mle. fit is just syntactic sugar i.e., we just thinly wrap fit_mle alongg with any initial steps one could want (e.g., maybe logic to init params or something). So in theory a user could also define fit for their specific IPP as well

[JoseKling]

So, I think this params_to_intensity and intensity_to_params is a bit weird. I guess we have two classes of intensities, parametric and non-parametric. You are implementing the fit method only for parametric intensity functions, so it might make sense to distinguish them. Just on the top of my head, not sure if this would be the bet way, but what if we do something like

abstract type ParametricIntensity end

struct ExponentialIntensity{R<:Real} <: ParametricIntensity
    a::R    # Could also store parameters as Vector or StaticVector and drop the 'mutable'
    b::R
end

StatsAPI.params(f::ExponentialIntensity) = (f.a, f.b)

Then in the negative_log_likelihood_ipp method we could replace the argument params_to_intensity with the F::Type{<:ParametericIntensity} and replace the first line where you build the intensity function with f = F(params...).

What do you think about this approach?

[rsenne]

So, I think this params_to_intensity and intensity_to_params is a bit weird

Yeah I agree something cleaner would be better. The one thing we will need to handle is how exactly we handle the parameters. The initial reason for my approach, was to handle how certain parameters get optimized. E.g., if a param needed to be postive we could optimize the log of the parameter. But we could do this maybe:

"""
Optimization parameters for ExponentialIntensity.

Optimize over θ = [log(a), b] so that a = exp(θ[1]) > 0.
"""
function StatsAPI.params(f::ExponentialIntensity)
    return [log(f.a), f.b]
end

"""
Construct ExponentialIntensity from optimization parameters θ.

θ[1] = log(a), θ[2] = b.
"""
function StatsAPI.from_params(::Type{ExponentialIntensity{R}}, θ) where {R<:Real}
    @assert length(θ) == 2
    a = exp(θ[1])
    b = θ[2]
    return ExponentialIntensity(R(a), R(b))
end
 """
    negative_loglikelihood_ipp(θ, ::Type{F}, times, tmin, tmax)

Negative log-likelihood for an inhomogeneous Poisson process with
parametric intensity `F` and event times `times` on [tmin, tmax].
"""
function negative_loglikelihood_ipp(θ, ::Type{F}, times, tmin, tmax) where {F<:ParametricIntensity}
    # Build intensity from optimization params
    f = StatsAPI.from_params(F, θ)

    # ∫ λ(t) dt over [tmin, tmax]
    Λ = integrated_intensity(f, tmin, tmax)

    # Sum of log intensities at event times
    log_sum = zero(eltype(θ))
    for t in times
        λ = f(t)
        if λ <= 0
            return typeof(log_sum)(Inf)  # invalid params
        end
        log_sum += log(λ)
    end

    return Λ - log_sum
end

[JoseKling]

Hmmm, now I understand the problem. The issue is when some prior transformation need to be done prior to the optimization.

I guess your idea works. We just define a default from_params method that should be overloaded only if needed. Perhaps calling it something like from_optimized_params would be clearer.

function from_optimized_params(PI::Type{<:ParametricIntensity}, θ)
    return PI(θ...)
end

[rsenne]

Okay so I tried to integrate our discussions a bit. In the most recent commit we do the following:

1. There is a new ParametricIntensity type.
2. To create a new ParametricIntensity one must specify a from_params and an initial_params function that will be called to transform params to the right domain and to give the optimizer a reasonable starting point. I'm not married to the names so please suggest what you like (you had suggested from_optimized_params but transform_params may also be okay.)
3. There is now an IntegralConfig so that users can specify whatever solvers/tolerances they desire for their specific use case.

Let me know which things aren't resolved/you'd like other changes on!

[JoseKling]

So, I was playing around a bit to see if I could find a way to simplify everything. I came up with what is the new branch pr-65-edits. Have a look and see if it makes sense or if it can be at least partly useful.

The main idea was to simplify the methods related to fit. What allowed this is that now you can pass the keyword arguments of the from_params method directly to fit. Also, the initial parameters must be provided by the user (not sure if this is ideal, but it does simplify things for PolynomialIntensity, which has variable number of parameters). So now we can do something like

h = History(rand(10), 0, 1, rand(10))

intensity = fit(PolynomialIntensity{Float64}, h, [1.0, 1.0, 1.0, 1.0]; link=:log)

ipp = fit(InhomogeneousPoissonProcess{PolynomialIntensity{Float64}, Uniform}, h, [30.0, 20.0]; link=:identity)

or

struct Constant{R} <: ParametricIntensity
a::R
end where {R<:Real}

(f::Constant, t) = f.a

# Might error because optimizer tests negative parameters, but works
fit(Constant, h, [4.0])

from_params(::Type{<:Constant}, params) = Constant(exp(params[1])

fit(Constant, h, [log(4.0)]) # No errors

There is actually no need to restrict fit to ParametricIntensity types, as long as the from_params method is implemented.

And you will see that I got rid of all the specific fit methods, but everything seems to work. Check if this is really the case.

[rsenne]

Okay, i looked at the branch and overall I agree the fit architecture has been improved. There are only a few remaining pinch points i think.

1.) PiecewiseConstant will need its own fit method. Given the function is piecewise--I'm not sure how one would even fit this using gradient descent. I would be open to either removing it entirely or retaining the separate fit method.

2.) I think if we want to allow things like identity links for methods like PolynomialIntensity we will need to keep the barrier method I previously had in the loss function, or alternatively we just only allow link functions like exp and softplus which will naturally keep the intensity in the right domain. If someone is fitting a polynomial anyway, my assumption is that they just want a smooth intensity function. In which case exp(lambda(t)) where lambda(t) is a polynomial shouldn't be a huge inconvenience. Or at the very least we can just make an informative error that gets thrown if we venture into such territory.

Thoughts?

[JoseKling]

I was not very thorough, so I apparently deleted/changed stuff I should not have.

1. In theory, if the breakpoints are fixed we could. But the closed-form log likelihood you implemented is, of course, much more efficient, so I would keep it.

2. We can keep the barriers you had before if it works. The thing is just that I am not completely sure it solves the problem (although it might avoid it most of the time). Let's say $f(\theta)$ is the objective function, then you are building a new objective as

$$g(\theta) = f(\theta) \text{, if } \lambda(t_i) \geq 0 \text{ for all } t_i$$ $$g(\theta) = \text{Inf, otherwise}$$

The thing is that the derivatives of $g$ are not defined in the region $\{\theta | \lambda(\theta, t_i) \leq 0 \text{ for any } t_i\}$.
I am not saying we should not do it, but I would think there are edge cases where this will present problems.

[rsenne]

Okay i re-added the PiecewiseConstant fit function and all tests pass again with your new interface for fitting. As for point 2, I agree that this isn't a perfect solution, but I think the alternative of not having to would be more prone to errors given the optimizer could wander into spaces where the intensity function has negative values. I think if our goal is to reasonably avoid all errors, then i think the solution is to make sure all Parametric models are fit as:

$$Poisson(\lambda_t(t) = \exp^{f(\theta)})$$

and not

$$Poisson(\lambda_t(t) = f(\theta))$$

I myself actually prefer this as if a user has a reason to fit some other model it leaves these concerns up to them.

[JoseKling]

I think if our goal is to reasonably avoid all errors, then i think the solution is to make sure all Parametric models are fit as:
$$Poisson(\lambda t(t) = expf(\theta))$$

I think this is reasonable, but I also think this is not an obvious behavior an user would expect. If I created an intensity, I would not expect it to be transformed in any way without me choosing it.

An alternative would be to set up something similar to the IntegrationConfig(). We could have a OptimizationConfig() for the to choose the solver, the constrains, or even add gradients. Not sure how exactly this would look like.

Anyway, I am fine leaving it as is or with this transformation you proposed for now. What you did is already a lot of functionality added, and we can improve this in the future.

[Copilot]

Pull request overview

This pull request adds comprehensive support for inhomogeneous Poisson processes (IPPs) to the point process modeling library. The implementation introduces time-varying intensity functions with both parametric and custom forms, numerical integration infrastructure, and maximum likelihood estimation with automatic differentiation for model fitting.

Key changes include:

• Five new parametric intensity function types (Polynomial, Exponential, Sinusoidal, PiecewiseConstant, LinearCovariate) with analytical optimizations where possible
• A flexible InhomogeneousPoissonProcess type that supports arbitrary callable intensity functions
• MLE-based fitting infrastructure using Optim.jl with automatic differentiation support via ForwardDiff

Reviewed changes

Copilot reviewed 11 out of 12 changed files in this pull request and generated 13 comments.

Show a summary per file

File	Description
Project.toml	Adds Integrals.jl and Optim.jl dependencies for numerical integration and optimization
test/Project.toml	Adds Optim.jl to test dependencies
src/PointProcesses.jl	Exports new IPP types, intensity functions, and configuration; includes new IPP modules
src/poisson/inhomogeneous/integration_config.jl	Defines IntegrationConfig for configurable numerical integration
src/poisson/inhomogeneous/parametric_intensity.jl	Trait interface for parametric intensity functions supporting MLE fitting
src/poisson/inhomogeneous/intensity_functions.jl	Implements five parametric intensity types and their parameter transformations
src/poisson/inhomogeneous/intensity_methods.jl	Provides analytical and numerical methods for integration and bounds computation
src/poisson/inhomogeneous/inhomogeneous_poisson_process.jl	Main IPP type implementing AbstractPointProcess interface with time-varying intensity
src/poisson/inhomogeneous/fit.jl	MLE fitting methods for parametric intensities and specialized fitting for piecewise constant
test/runtests.jl	Adds test suite inclusion for inhomogeneous Poisson processes
test/inhomogeneous_poisson_process.jl	Comprehensive tests covering all intensity types, fitting, simulation, and interface methods
docs/src/api.md	Documents new IPP types, intensity functions, and configuration in API reference

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[rsenne]

Okay all tests now pass. Let me know if there is anything else you would like

[JoseKling]

Only problem I see are the intensity_bound methods. Left comments on each of the methods.

The rest looks fine.

[rsenne]

Should be all done @JoseKling let me know if anything else pops up

[JoseKling]

That's good. I just suggested not having the min_lookahead keyword. I just don't see why we should have it. If the length of the process is of this order, I don't see the problem of letting lookahead be smaller than 1e-3, especially since the user cannot control this in the simulate_ogata method.

I am approving the PR, just either accept the suggestions or, if I am overlooking something, ignore them.