X-DataInitiative / tick

Module for statistical learning, with a particular emphasis on time-dependent modelling
https://x-datainitiative.github.io/tick/
BSD 3-Clause "New" or "Revised" License
486 stars 105 forks source link

Negative amplitudes in `HawkesSumGaussians` #467

Open trouleau opened 3 years ago

trouleau commented 3 years ago

Hello,

I believe there is a bug in the implementation of the algorithm HawkesSumGaussians. In the paper defining the algorithm, the amplitudes of the kernel bases are defined to be non-negative (See Equation 4). However, the amplitudes attributes of a fitted HawkesSumGaussians can be negative.

Here is a minimal example to reproduce the issue (tested on the current version tick==0.7.0.1):

import numpy as np
import tick.hawkes

# Generate model parameters
n_nodes = 4
adjacency = np.array([[1, 1, 0, 0],
                      [0, 1, 0, 0],
                      [0, 0, 1, 1],
                      [0, 0, 0, 1]]) * 0.5
baseline = np.ones(n_nodes) * 0.1
decays = 1.0

# Simulate a realization of the process
hawkes_simu = tick.hawkes.SimuHawkesExpKernels(adjacency=adjacency, decays=decays,
                                               baseline=baseline, max_jumps=50e3,
                                               verbose=False, seed=123)
hawkes_simu.simulate()

# Run inference with MLE-SGLP
hawkes = tick.hawkes.HawkesSumGaussians(max_mean_gaussian=10.0, n_gaussians=5,)
hawkes.fit(hawkes_simu.timestamps)

print(hawkes.hawkes.amplitudes.min())

I believe the issue comes from the definition of the soft-thresholding function soft_thres. In the code, the sign function is defined on {-1, 1}, but I think it should be defined on {0, 1}, i.e., it should be zero if the input is negative (not -1).

If you agree that this is indeed a bug, then I guess the fix should be as easy as replacing the soft_thres function into

double soft_thres(double z, double alpha) {
  return (z > 0 ? std::max(std::abs(z) - alpha, 0.) : 0.);
}
Jbollenbacher commented 3 years ago

I have noticed this problem too.

Edit: Having looked into this a bit more, there seems to be no lower limit on amplitude, and you can end up with very large negative amplitudes in the estimated intensity functions. This isn't reasonable mathematically, so I would say this is definitely a bug.