Hakes processes - random time change and QQ plots (#497)

Hawkes processes - random time change and QQ plots

.gitignore

@@ -56,4 +56,5 @@ tools/benchmark/data
 
 # Virtual Environments
 env*/
-venv*/
+venv*/
+tickf.yml

examples/plot_hawkes_basis_kernels.py

@@ -38,8 +38,13 @@ def g2(t):
     return np.cos(np.pi * (t / 10 + 1)) + 1.1
 
 
+############################################################################
+# Simulate
+############################################################################
 t_values = np.linspace(0, 20, 1000)
-u_values = [(0.007061, 0.001711), (0.005445, 0.003645), (0.003645, 0.005445),
+u_values = [(0.007061, 0.001711),
+            (0.005445, 0.003645),
+            (0.003645, 0.005445),
             (0.001790, 0.007390)]
 
 hawkes = SimuHawkes(baseline=[1e-5, 1e-5], seed=1093, verbose=False)
@@ -51,17 +56,29 @@ def g2(t):
 
 hawkes.end_time = end_time
 hawkes.simulate()
-ticks = hawkes.timestamps
 
-# And then perform estimation with two basis kernels
+
+#############################################################################
+# Fit
+#############################################################################
+timestamps = hawkes.timestamps
 kernel_support = 20
 n_basis = 2
+end_time = 1e9
+C = 1e-3
+kernel_size = 40
+max_iter = 100
 
 em = HawkesBasisKernels(kernel_support, n_basis=n_basis,
                         kernel_size=kernel_size, C=C, n_threads=4,
                         max_iter=max_iter, verbose=False, ode_tol=1e-5)
-em.fit(ticks)
+em.fit(timestamps)
+
 
+#############################################################################
+# Plot
+#############################################################################
+show = True
 fig = plot_hawkes_kernels(em, hawkes=hawkes, support=19.9, show=False)
 for ax in fig.axes:
     ax.set_ylim([0, 0.025])
@@ -70,4 +87,5 @@ def g2(t):
 for ax in fig.axes:
     ax.set_ylim([0, 0.5])
 
-plt.show()
+if show:
+    plt.show()

examples/plot_hawkes_estimated_intensity.py

@@ -1,7 +1,8 @@
 """
-============================================
-Plot estimated intensity of Hawkes processes
-============================================
+================================================
+Plot estimated intensity of Hawkes processes and
+assess goodness of fit via QQ plots
+=================================================
 
 This examples shows how the estimated intensity of a learned Hawkes process
 can be plotted. In this example, the data has been generated so we are able
@@ -11,32 +12,60 @@
 
 import matplotlib.pyplot as plt
 
-from tick.hawkes import SimuHawkesSumExpKernels, HawkesSumExpKern
-from tick.plot import plot_point_process
+from tick.hawkes import (
+    SimuHawkesSumExpKernels,
+    HawkesSumExpKern
+)
+from tick.hawkes import SimuHawkesExpKernels  # NOQA
+from tick.hawkes import HawkesExpKern  # NOQA
+from tick.plot import plot_point_process, qq_plots
 
-end_time = 1000
 
+##########################################################################
+# simulate
+##########################################################################
+Simulator = SimuHawkesSumExpKernels
+end_time = 1000
 decays = [0.1, 0.5, 1.]
 baseline = [0.12, 0.07]
 adjacency = [[[0, .1, .4], [.2, 0., .2]], [[0, 0, 0], [.6, .3, 0]]]
+model = Simulator(
+    adjacency=adjacency,
+    decays=decays,
+    baseline=baseline,
+    end_time=end_time,
+    verbose=False,
+    seed=1039,
+)
+model.track_intensity(0.1)
+model.simulate()
 
-hawkes_exp_kernels = SimuHawkesSumExpKernels(
-    adjacency=adjacency, decays=decays, baseline=baseline, end_time=end_time,
-    verbose=False, seed=1039)
 
-hawkes_exp_kernels.track_intensity(0.1)
-hawkes_exp_kernels.simulate()
+##########################################################################
+# fit
+##########################################################################
+timestamps = model.timestamps
+Fitter = HawkesSumExpKern
+decays = [0.1, 0.5, 1.]
+kwargs = {}
+if Fitter == HawkesSumExpKern:
+    if 'penalty' not in kwargs:
+        kwargs['penalty'] = 'elasticnet'
+        kwargs['elastic_net_ratio'] = 0.8
+learner = Fitter(decays=decays, **kwargs)
+learner.fit(timestamps)
 
-learner = HawkesSumExpKern(decays, penalty='elasticnet', elastic_net_ratio=0.8)
-learner.fit(hawkes_exp_kernels.timestamps)
 
+##########################################################################
+# plot intensities
+##########################################################################
 t_min = 100
 t_max = 200
+show = True
 fig, ax_list = plt.subplots(2, 1, figsize=(10, 6))
-learner.plot_estimated_intensity(hawkes_exp_kernels.timestamps, t_min=t_min,
+learner.plot_estimated_intensity(model.timestamps, t_min=t_min,
                                  t_max=t_max, ax=ax_list)
-
-plot_point_process(hawkes_exp_kernels, plot_intensity=True, t_min=t_min,
+plot_point_process(model, plot_intensity=True, t_min=t_min,
                    t_max=t_max, ax=ax_list)
 
 # Enhance plot
@@ -54,5 +83,39 @@
 
     ax.legend()
 
-fig.tight_layout()
-plt.show()
+if show:
+    fig.tight_layout()
+    plt.show()
+
+
+def simulated_v_estimated_qq_plots(
+        model,
+        learner,
+        show=True,
+):
+    fig, ax_list = plt.subplots(2, 1, figsize=(10, 6))
+    timestamps = model.timestamps
+    end_time = model.end_time
+    learner.qq_plots(events=timestamps, end_time=end_time, ax=ax_list)
+    model.store_compensator_values()
+    qq_plots(model, ax=ax_list)
+
+    # Enhance plot
+    for ax in ax_list:
+        # Set labels to both plots
+        ax.lines[0].set_label('estimated')
+        ax.lines[2].set_label('original')
+
+        # Change original intensity style
+        ax.lines[2].set_alpha(0.6)
+        ax.lines[3].set_alpha(0.6)
+
+        ax.lines[2].set_markerfacecolor('orange')
+        ax.lines[2].set_markeredgecolor('orange')
+
+        ax.legend()
+
+    if show:
+        fig.tight_layout()
+        plt.show()
+    return fig

examples/plot_hawkes_time_func_simu.py

@@ -11,8 +11,11 @@
 
 from tick.base import TimeFunction
 from tick.hawkes import SimuHawkes, HawkesKernelExp, HawkesKernelTimeFunc
-from tick.plot import plot_point_process
+from tick.plot import plot_point_process, qq_plots as _qq_plots
 
+###############################################################################
+# instantiate
+###############################################################################
 t_values = np.array([0, 1, 1.5], dtype=float)
 y_values = np.array([0, .2, 0], dtype=float)
 tf1 = TimeFunction([t_values, y_values],
@@ -25,17 +28,26 @@
                    dt=0.1)
 kernel_2 = HawkesKernelTimeFunc(tf2)
 
-hawkes = SimuHawkes(
+model = SimuHawkes(
     kernels=[[kernel_1, kernel_1], [HawkesKernelExp(.07, 4), kernel_2]],
     baseline=[1.5, 1.5], verbose=False, seed=23983)
 
+
+###############################################################################
+# simulate
+###############################################################################
 run_time = 40
 dt = 0.01
-hawkes.track_intensity(dt)
-hawkes.end_time = run_time
-hawkes.simulate()
-
-fig, ax = plt.subplots(hawkes.n_nodes, 1, figsize=(14, 8))
-plot_point_process(hawkes, t_max=20, ax=ax)
-
+model.track_intensity(dt)
+model.end_time = run_time
+model.simulate()
+
+
+###############################################################################
+# plot
+###############################################################################
+fig1, ax = plt.subplots(model.n_nodes, 1, figsize=(14, 8))
+plot_point_process(model, t_max=20, ax=ax)
+model.store_compensator_values()
+fig2 = _qq_plots(model, show=False)
 plt.show()

examples/plot_hawkes_varying_baseline.py

@@ -13,11 +13,14 @@
 """
 
 import matplotlib.pyplot as plt
-
-from tick.plot import plot_hawkes_baseline_and_kernels
+from tick.plot import plot_hawkes_baseline_and_kernels, qq_plots
 from tick.hawkes import (SimuHawkesSumExpKernels, SimuHawkesMulti,
                          HawkesSumExpKern)
 
+
+################################################################################
+# instantiate and run
+################################################################################
 period_length = 300
 baselines = [[0.3, 0.5, 0.6, 0.4, 0.2, 0], [0.8, 0.5, 0.2, 0.3, 0.3, 0.4]]
 n_baselines = len(baselines[0])
@@ -26,22 +29,61 @@
 
 # simulation
 hawkes = SimuHawkesSumExpKernels(baseline=baselines,
-                                 period_length=period_length, decays=decays,
-                                 adjacency=adjacency, seed=2093, verbose=False)
+                                 period_length=period_length,
+                                 decays=decays,
+                                 adjacency=adjacency,
+                                 seed=2093,
+                                 verbose=False,
+                                 )
 hawkes.end_time = 1000
 hawkes.adjust_spectral_radius(0.5)
-
 multi = SimuHawkesMulti(hawkes, n_simulations=4)
 multi.simulate()
 
 # estimation
 learner = HawkesSumExpKern(decays=decays, n_baselines=n_baselines,
                            period_length=period_length)
-
 learner.fit(multi.timestamps)
 
+
+################################################################################
 # plot
+################################################################################
+show = True
 fig = plot_hawkes_baseline_and_kernels(learner, hawkes=hawkes, show=False)
-fig.tight_layout()
+if show:
+    fig.tight_layout()
+    plt.show()
+
+
+def simulated_v_estimated_qq_plots(
+        model,
+        learner,
+        show=True,
+):
+    fig, ax_list = plt.subplots(2, 1, figsize=(10, 6))
+    timestamps = model.timestamps
+    end_time = model.end_time
+    learner.qq_plots(events=timestamps, end_time=end_time, ax=ax_list)
+    model.store_compensator_values()
+    qq_plots(model, ax=ax_list)
+
+    # Enhance plot
+    for ax in ax_list:
+        # Set labels to both plots
+        ax.lines[0].set_label('estimated')
+        ax.lines[2].set_label('original')
+
+        # Change original intensity style
+        ax.lines[2].set_alpha(0.6)
+        ax.lines[3].set_alpha(0.6)
+
+        ax.lines[2].set_markerfacecolor('orange')
+        ax.lines[2].set_markeredgecolor('orange')
+
+        ax.legend()
 
-plt.show()
+    if show:
+        fig.tight_layout()
+        plt.show()
+    return fig

examples/plot_hawkes_varying_baseline_simulation.py

@@ -11,8 +11,12 @@
 from tick.base import TimeFunction
 from tick.hawkes import SimuHawkesExpKernels
 
-from tick.plot import plot_point_process
+from tick.plot import plot_point_process, qq_plots as _qq_plot
 
+
+###############################################################################
+# instantiate
+##############################################################################
 period_length = 100
 t_values = np.linspace(0, period_length)
 y_values = 0.2 * np.maximum(
@@ -22,24 +26,49 @@
 
 decay = 0.1
 adjacency = np.array([[0.5]])
+hawkes = SimuHawkesExpKernels(adjacency,
+                              decay,
+                              baseline=baselines,
+                              seed=2093,
+                              verbose=False,
+                              )
+
 
-hawkes = SimuHawkesExpKernels(adjacency, decay, baseline=baselines, seed=2093,
-                              verbose=False)
+###############################################################################
+# simulate
+##############################################################################
+period_length = 100
 hawkes.track_intensity(0.1)
 hawkes.end_time = 6 * period_length
 hawkes.simulate()
 
-fig, ax = plt.subplots(1, 1, figsize=(10, 4))
 
+###############################################################################
+# plot intensity
+##############################################################################
+fig1, ax = plt.subplots(1, 1, figsize=(10, 4))
 plot_point_process(hawkes, ax=ax)
-
 t_values = np.linspace(0, hawkes.end_time, 1000)
-ax.plot(t_values, hawkes.get_baseline_values(0, t_values), label='baseline',
+ax.plot(t_values,
+        hawkes.get_baseline_values(0, t_values),
+        label='baseline',
         ls='--', lw=1)
 ax.set_ylabel("$\lambda(t)$", fontsize=18)
 ax.legend()
-
 plt.title("Intensity Hawkes process with exponential kernel and varying "
           "baseline")
-fig.tight_layout()
+
+
+###############################################################################
+# qq plot
+##############################################################################
+hawkes.store_compensator_values()
+fig2 = _qq_plot(hawkes)
+
+
+###############################################################################
+# show
+##############################################################################
+fig1.tight_layout()
+fig2.tight_layout()
 plt.show()