From 0d0bdb305f3e873f70f52bf2fd48aa1da24f2739 Mon Sep 17 00:00:00 2001 From: Levi John Wolf Date: Mon, 13 Oct 2025 16:01:20 +0100 Subject: [PATCH 1/2] add start of hawkes MLE --- pointpats/_hawkes.py | 76 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 76 insertions(+) create mode 100644 pointpats/_hawkes.py diff --git a/pointpats/_hawkes.py b/pointpats/_hawkes.py new file mode 100644 index 00000000..7349cc6d --- /dev/null +++ b/pointpats/_hawkes.py @@ -0,0 +1,76 @@ +# Estimate a spatio-temporal Hawkes processes using maximum likelihood estimation +import numpy +from scipy.spatial import distance + +class UnivariateStationaryHawkesProcess(): + def __init__(self, geo_kernel, time_kernel): + self.geo_kernel = spatial_kernel, + self.time_kernel = time_kernel + + + def _triggering_kernel_st( + s: numpy.ndarray, + t: float, + w: numpy.ndarray, + u: float, + baseline_excitation = 0, + time_lengthscale = 1, + space_lengthscale = 1 + ) -> float: + if t < u: + return 0.0 + return ( + numpy.log(baseline_excitation) + ( + # temporal, should be the log of the kernel value + lumpy.log(self.time_kernel(t, u, time_lengthscale)) + ) + ( + # spatial, should be the log of the kernel value + numpy.log(self.spatial_kernel(s, w, space_lengthscale)) + ) + ) + + def _triggering_kernel_delta( + self, + distance: numpy.ndarray, + duration: numpy.ndarray, + baseline_excitation = 0, + time_lengthscale = 1, + space_lengthscale = 1, + ): + if duration <= 0: + return 0 +` return numpy.log(baseline_excitation) + ( + # temporal + - numpy.log(time_lengthscale) - duration/time_lengthscale + ) + ( + # spatial + -numpy.log(2*numpy.pi*space_lengthscale**2) - distance/(2*space_lengthscale**2) + ) + + def _log_likelihood( + baseline_excitation, + time_lengthscale, + space_lengthscale, + mu, + mu0, + distances, + durations + ): + # triggering function component + n_samples,_ = distances.shape + for i in range(1, n) + for j in range(i, n): + if durations[i,j] > 0: + ll[i] += self._triggering_kernel_delta( + distances[i,j], + durations[i,j], + baseline_excitation=baseline_excitation, + time_lengthscale=time_lengthscale, + space_lengthscale=space_lengthscale + ) + # integral component + + + def fit(self, geometry, times): + # maximize log-likelihood over baseline_excitation, time_lengthscale, space_lengthscale + \ No newline at end of file From 28ab8dbef4433340c4e87e492e2b2adc6d8a8124 Mon Sep 17 00:00:00 2001 From: Levi John Wolf Date: Mon, 13 Oct 2025 16:55:51 +0100 Subject: [PATCH 2/2] add multivariate hawkes process stub --- pointpats/_hawkes.py | 114 +++++++++++++++++++++++++++++++++++-------- 1 file changed, 95 insertions(+), 19 deletions(-) diff --git a/pointpats/_hawkes.py b/pointpats/_hawkes.py index 7349cc6d..dc5ef3bc 100644 --- a/pointpats/_hawkes.py +++ b/pointpats/_hawkes.py @@ -1,12 +1,16 @@ # Estimate a spatio-temporal Hawkes processes using maximum likelihood estimation import numpy from scipy.spatial import distance +from scipy import optimize -class UnivariateStationaryHawkesProcess(): - def __init__(self, geo_kernel, time_kernel): - self.geo_kernel = spatial_kernel, - self.time_kernel = time_kernel +# maybe the better way to do this would be to implement +# the stan code, then allow the user to select various kernels? +# otherwise we'll have a lot of very long classes. +class UnivariateStationarySeparableSpatialHawkes(): + def __init__(self, space_kernel, time_kernel): + self.space_kernel = spatial_kernel, + self.time_kernel = time_kernel def _triggering_kernel_st( s: numpy.ndarray, @@ -17,7 +21,7 @@ def _triggering_kernel_st( time_lengthscale = 1, space_lengthscale = 1 ) -> float: - if t < u: + if t <= u: return 0.0 return ( numpy.log(baseline_excitation) + ( @@ -25,7 +29,7 @@ def _triggering_kernel_st( lumpy.log(self.time_kernel(t, u, time_lengthscale)) ) + ( # spatial, should be the log of the kernel value - numpy.log(self.spatial_kernel(s, w, space_lengthscale)) + numpy.log(self.space_kernel(s, w, space_lengthscale)) ) ) @@ -36,32 +40,33 @@ def _triggering_kernel_delta( baseline_excitation = 0, time_lengthscale = 1, space_lengthscale = 1, - ): + ) -> float: if duration <= 0: - return 0 -` return numpy.log(baseline_excitation) + ( + return 0.0 + return numpy.log(baseline_excitation) + ( # temporal - - numpy.log(time_lengthscale) - duration/time_lengthscale + numpy.log(self.time_kernel(duration, time_lengthscale)) ) + ( # spatial - -numpy.log(2*numpy.pi*space_lengthscale**2) - distance/(2*space_lengthscale**2) + numpy.log(self.space_kernel(distance, space_lengthscale)) ) def _log_likelihood( baseline_excitation, time_lengthscale, space_lengthscale, - mu, - mu0, + intensity, + baseline_intensity, distances, durations - ): + ) -> float: # triggering function component n_samples,_ = distances.shape + excitation = numpy.zeros(n_samples) for i in range(1, n) for j in range(i, n): if durations[i,j] > 0: - ll[i] += self._triggering_kernel_delta( + excitation[i] += self._triggering_kernel_delta( distances[i,j], durations[i,j], baseline_excitation=baseline_excitation, @@ -69,8 +74,79 @@ def _log_likelihood( space_lengthscale=space_lengthscale ) # integral component - + baseline_all = baseline_intensity * self.space_window * self.time_window + + baseline_time_nospace = baseline_excitation * self.time_kernel(self.waits_, time_lengthscale).sum() + return -( + excitation.sum() + + self.integral_value*intensity + + baseline_all + + baseline_time_nospace - + self.integral_value + ) - def fit(self, geometry, times): - # maximize log-likelihood over baseline_excitation, time_lengthscale, space_lengthscale - \ No newline at end of file + def fit(self, geometry: geopandas.GeoSeries, times: pandas.Series[pd.Timestamp]) -> UnivariateStationaryHawkesProcess: + coords = geometry.get_coordinates() + durations = numpy.subtract.outer(times, times) + distances = distance.cdist(coords, coords) + # we should probably estimate on the unit square+hour and then re-transform? + self.time_window_ = durations.max() - durations.min() + + self.space_window_ = distances.max() - distances.min() + self.waits_ = self.time_window_ - times # now time is denominated as float timestep since start at 0 + # maximize log-likelihood over mu, mu0, baseline_excitation, time_lengthscale, space_lengthscale + self._calcluate_integral() # calculate the integral value over space and time + + def score(vars): + return - self._log_likelihood( + vars[0], # baseline_excitation + vars[1], # time_lengthscale + vars[2], # space_lengthscale + vars[3], # intensity + vars[4], # baseline_intensity + distances, + durations + ) + + result = optimize.minimize( + score, + x0 = numpy.array([0.1, 1.0, 1.0, 0.1, 0.1]), # initial guess + args = ( + intensity, + baseline_intensity, + distances, + durations + ), + bounds = ( + (1e-5, None), # baseline excitation + (1e-5, None), # time lengthscale + (1e-5, None), # space lengthscale + (1e-5, None), # intensity + (1e-5, None) # baseline intensity + ), + method = 'L-BFGS-B' + ) + self.baseline_excitation_ = result.x[0] + self.time_lengthscale_ = result.x[1] + self.space_lengthscale_ = result.x[2] + self.intensity_ = result.x[3] + self.baseline_intensity_ = result.x[4] + return self + + def _calculate_integral(self): + # calculate the integral of the triggering function over space and time + # I think if proper kernels are used, this integral is n? + self.integral_value = len(self.waits_) # placeholder, should be the integral value only if proper kernels are used + +class MultivariateStationarySeparableSpatialHawkes(): + def __init__(self, space_kernel, time_kernel): + self.space_kernel = spatial_kernel, + self.time_kernel = time_kernel + + def fit(self, geometry: geopandas.GeoSeries, times: pandas.Series[pd.Timestamp], types: pandas.Series[int]) -> MultivariateStationaryHawkesProcess: + # similar to univariate but now we have a matrix of baseline excitations and intensities. + # Basically, we do the same as in the univariate, but the triggering loop has to be done + # across all pairs of event types. Think of the univariate case above as a within-type triggering equation, + # while doing this from type A (n_samples_A) to type B (n_samples_B) would be + # a between-type triggering equation, iterating for i in range(n_samples_A) and j in range(n_samples_B), calculating + # using the cross-type matrix of distances/durations + pass \ No newline at end of file