Add parametric baseline intensity?

[kochbj]

Hi Dean,

This looks like an awesome package and I'm excited to try it out! I'm interested in modeling citations with a Hawkes process where the baseline rate decays according to an exponential kernel as well, Something like:

$\mu_t= \alpha * \exp(-t * \delta) $

Any advice on how I might adapt the code to do something like this? Perhaps without estimating the delta parameter since I assume that complicates things?

Thanks in advance,
Bernie

[dm13450]

In short, yes, what you've described is exactly what this package can do with a little bit of work. It's been on my to-do list to write about how to do it.

If alpha and delta are known and fixed, then it should work out of the box. I haven't tested it, but this code will sample the kappa and kernel parameters.

function fit(events::Array{<:Number, 1}, maxT::Number, its::Int64)

    nEvents = length(events)
    eventDifferences = event_difference_list(events)
    maxTDifferences = maxT .- events

    kappaSamples = zeros(Float64, its)
    kernSamples = zeros(Float64, its)

    bgSamples[1] = rand()
    kappaSamples[1] = rand()
    kernSamples[1] = rand()

    kernSample = kernSamples[1]
    kernDistribution::Exponential{Float64} = Distributions.Exponential(1/kernSample)
    kernFunc(x::Number) = Distributions.pdf(kernDistribution, x)

   bgFunc(t) = 0.01 * exp(-0.01 * t) 

    parentVectorSample, numBG, chEvents, shiftTimes = sample_parents(events, bgFunc, kappaSamples[1], kernFunc, eventDifferences)

    for i = 2:its

        H_tilde::Float64 = sum(cdf.(kernDistribution, maxTDifferences))
        kappaSamples[i] = Distributions.rand(Distributions.Gamma(0.01 + sum(chEvents), 1/(0.01 + H_tilde)))

        kernSample = Distributions.rand(Distributions.Gamma(0.01 + length(shiftTimes), 1/(0.01 + sum(shiftTimes))))
        kernSamples[i] = kernSample
        kernDistribution = Distributions.Exponential(1/kernSamples[i])

        parentVectorSample, numBG, chEvents, shiftTimes = sample_parents(events, bgFunc, kappaSamples[i], kernFunc, eventDifferences)
    end
    kappaSamples, kernSamples
end

To estimate alpha and delta you'll need to write a sampler that takes in the background events and returns a sample from the alpha and delta posterior.

function sample_new_alpha_delta(bgEvents)

  # do some sampling

  return alpha, delta
end

you might find http://dm13450.github.io/2020/11/03/BayesPointProcess.html this helpful for how to estimate the background rate.

Once you've got that function, its a case of including it in the fit function such that with each iteration you create a new bgFunc with the latest sample of the parameters.

Let me know if you need anymore help!

[kochbj]

Thanks so much for the fast and thorough response. Really appreciate it! Since I have very little experience with either point processes or Julia, I was wondering if you had any further advice on the best way to sample. Is there a Gibbs sampling approach for this parameterization that makes sense? If not, would you just write a simple MCMC sampler into this or is there a more faster/cleaner implementation in Julia that I could just plug in?

Best,
B

[dm13450]

Ok so I've created this repository that contains a notebook that implements the basic version of your background rate:
https://github.com/dm13450/Hawkes-Processes-Examples/blob/main/Exponential%20Decay%20in%20Background%20.ipynb

It's a little rough and ready, but should act as a good springboard for you to work off.

Again, happy to help/answer any more questions you have on this type of work.

[kochbj]

Hi Dean,

Sorry for the slow response. This is totally great and definitely a sufficient springboard. Appreciate the spoonfeeding and hopefully the example will be helpful for somebody else. Will follow up if I have any more qs.

Thanks again,
B