Remarks about package [JuliaCon Proceedings review]

[gdalle]

You know the drill by now. Here are a few general questions I had while reading the paper and docs:

https://github.com/JuliaCon/proceedings-review/issues/133

Functionalities

• The main reason I wrote PointProcesses.jl was to work with processes that generate arbitrary marks (not just multivariate counting processes). I think it is possible to do that with JumpProcesses.jl but I'd like some pointers
• Other interesting aspects of my library include loglikelihood computation and fitting. I think fitting is not in your scope, but is there any chance you would consider implementing a loglikelihood computation given a process and a sample?

[isaacsas]

Couldn't mark sampling be encoded in the affect! function? You can always setup your parameter object to have a place to store the marks (if you don't want to store them in the solution u). We should probably add a doc example showing how to do this (and how to get at the rng we are using so user samplers can also use it).

We don't have any explicit higher level fitting functionality here, but we should compose with higher level packages in SciML and such for fitting (modulo that AD can't really be used through the aggregation methods -- but this is something we are working on). We could certainly consider adding such computations here if it makes sense. Do you have any specific references for the case you are interested in?

[ChrisRackauckas]

Couldn't mark sampling be encoded in the affect! function? You can always setup your parameter object to have a place to store the marks (if you don't want to store them in the solution u). We should probably add a doc example showing how to do this (and how to get at the rng we are using so user samplers can also use it).

Yes, the marks are just handled in the affect! functions. The only thing that needs specialization there are RegularJumps, and I believe there's stuff written on that aspect already?

Other interesting aspects of my library include loglikelihood computation and fitting. I think fitting is not in your scope, but is there any chance you would consider implementing a loglikelihood computation given a process and a sample?

The general loglikelihoods can only be defined through sampling. Indeed the simple cases can have analytical solution and it can make sense to write that down, but for as they become general there's no analytical solution and sampling 10,000 times is the way to estimate the likelihood.

[isaacsas]

Yes, RegularJumps should support marks with the newer interface (I think the person who added that didn't add any docs unfortunately).

[gzagatti]

Since a marked point process is just a ground point process with a mark distribution. It should be pretty straightforward to implement marked point processes by drawing from an arbitrary distribution inside of the affect! function.

Yes, the marks are just handled in the affect! functions. The only thing that needs specialization there are RegularJumps, and I believe there's stuff written on that aspect already?

As an example, I put together a marked TPP using RegularJump in the past. This uses the old interface, but it should give an idea of what's possible.

function spair_prob(g::Graphs.AbstractSimpleGraph, tspan, p)

    e = ne(g)
    ē = convert(eltype(e), nv(g) * (nv(g) - 1) / 2 - ne(g))

    λ₁, λ₂ = p
    λg = λ₁ * e + λ₂ * ē

    d = MixtureModel(
        DiscreteUniform[DiscreteUniform(1, e), DiscreteUniform(e + 1, e + ē)],
        [λ₁ * e / λg, λ₂ * ē / λg],
    )
    # a = vcat(fill(1, e), fill(0, ē))
    # p = a.*(λ₁/λg) .+ (1 .- a).*(λ₂/λg)
    # d = Categorical(p)

    u0 = zeros(e + ē)

    dprob = DiscreteProblem(u0, tspan, [λg])

    jump = RegularJump(spair_rate, spair_count, 1; mark_dist = d)
    jprob = JumpProblem(dprob, Direct(), jump)

    return jprob

end

The general loglikelihoods can only be defined through sampling. Indeed the simple cases can have analytical solution and it can make sense to write that down, but for as they become general there's no analytical solution and sampling 10,000 times is the way to estimate the likelihood.

I might have a different understanding of the loglikelihood. But as I see, we have that the loglikelihood of any TPP which satisfy some basic assumptions is as following.

$$ \begin{split} \ell &(H{T^-}) = \ &= \sum{n=1}^N \left( \log \lambda^\ast(t_n, k_n) \right) - \int0^T \sum{k \in \mathcal{K}} \lambda^\ast (u, k) du \ &= \sum{n=1}^N \left( \log \frac{ \partial \Lambda{k_n}^\ast}{ \partial t_n} (t_n) \right) - \Lambda_g^\ast (T) \end{split} $$

Where $\lambda$ represents the conditional intensity and $\Lambda$ the compensator. The JumpProcesses.jl refers to the conditional intensity as a rate.

Since all the aggregators in JumpProcesses.jl expect the user to supply the intensity function, it would be possible to compute the log-likelihood as following:

1. Simulate the path
2. Given the path, integrate the conditional intensity numerically to obtain the compensator. Or compute the integrator if the user supply a closed form.
3. Compute the log-likelihood with the formula above

Does it make sense?

The main reason I wrote PointProcesses.jl was to work with processes that generate arbitrary marks (not just multivariate counting processes). I think it is possible to do that with JumpProcesses.jl but I'd like some pointers

I gave a look at the PointProcesses.jl documentation but didn't find any tutorials or examples that I could port to JumpProcesses.jl. If you have something please share. Alternatively, I could draft a simple example that walks the user through a simple TPP problem.

@gdalle let me know if that would be enough to address the issues you raised in the original post.

[gdalle]

Hey there, The new advanced TPP tutorial is great and does answer my questions about marks! I didn't read it in full detail but from the looks of it, you don't really exploit the PointProcesses.jl interface (and that's okay). My package is mostly designed to

• get free loglikelihood and simulation from intensity functions.
• implement a few simple Poisson processes

Since you rebuild both loglikelihood and simulation from scratch, I think it's not worth trying to stick to my design, especially since it's old and unmaintained. A small mention of my package (as well as others in the Julia ecosystem) in the paper or docs would be more than enough. What do you think?

[gzagatti]

@gdalle thanks for your suggestion. I have opened a PR that resolves this issue. Let us know what you think.

In the tutorial, I defined a concrete implementation for all the methods that you defined in abstract_point_process.jl. I liked the approach that you create and I thought that it would fit well with a SciML implementation of TPPs. In particular, I liked your approach to handling the event history. A long time ago, I started a private TPP package which I used to explore basic TPP models and had the same flavour as yours but over-complicated the event history.

My first version of the tutorial got the loglikelihood for free from your package. But then I realized that it would be more efficient to defined a custom implementation. Given the history of events, SciMLPointProcess solves a DiscreteProblem from OrdinaryDiffEq which returns an object with all the stopping times of the intensity path. This takes O(N). Therefore, it would be inefficient to have another loop that computes the log_intensity for each stopping time as in the implementation of your package, because we would have to re-compute the path of the intensity for each call to log_intensity which would be O(N^2). Most of my additions to your API were to allow the methods to take a pre-computed intensity path. Otherwise, I used the same signatures.

I also modified the simulation part because this is exactly what JumpProcesses.jl is meant to do.

I hope that clarifies your issues.

[gdalle]

If you like the approach, by all means you can keep it, I just meant you shouldn't feel obliged :) And I'm glad you came up with a more effcient loglikelihood evaluation. That might actually inspire changes to my package if I ever revisit it.

I'll read the revamped tutorials and paper in the coming days