From 84456aab447a49e2d1c7bc96a5029d97993846a4 Mon Sep 17 00:00:00 2001
From: Ryan Senne <50930199+rsenne@users.noreply.github.com>
Date: Fri, 8 May 2026 22:19:38 -0400
Subject: [PATCH 1/6] Add time rescaling to inhomoggeneous

---
 docs/examples/Inhomogeneous.jl                | 20 ++++++
 .../inhomogeneous_poisson_process.jl          | 21 ++++++
 test/hypothesis_tests.jl                      | 37 ++++++++++
 test/inhomogeneous_poisson_process.jl         | 70 +++++++++++++++++++
 4 files changed, 148 insertions(+)

diff --git a/docs/examples/Inhomogeneous.jl b/docs/examples/Inhomogeneous.jl
index 256e763..c271e5b 100644
--- a/docs/examples/Inhomogeneous.jl
+++ b/docs/examples/Inhomogeneous.jl
@@ -215,6 +215,26 @@ scatter!(
     alpha=0.6,
 )
 
+# ## Goodness-of-Fit via Time Rescaling
+# By the time-rescaling theorem, if a fitted intensity ``\hat{λ}(t)`` is correct, then the
+# events transformed by the compensator ``Λ(t) = ∫_{t_{\min}}^{t} \hat{λ}(s) ds`` form a
+# unit-rate Poisson process. We can use `MonteCarloTest` with `KSDistance{Exponential}` to
+# check whether the transformed inter-event times look Exp(1):
+pp_poly_full = InhomogeneousPoissonProcess(pp_poly)
+pp_gauss_full = InhomogeneousPoissonProcess(pp_gauss)
+
+mc_piecewise = MonteCarloTest(KSDistance{Exponential}, pp_piecewise, h; n_sims=500, rng=rng)
+mc_poly = MonteCarloTest(KSDistance{Exponential}, pp_poly_full, h; n_sims=500, rng=rng)
+mc_gauss = MonteCarloTest(KSDistance{Exponential}, pp_gauss_full, h; n_sims=500, rng=rng)
+
+println("\nGoodness-of-fit p-values (Monte Carlo KS test):") # hide
+println("  Piecewise Constant: ", round(pvalue(mc_piecewise); digits=3)) # hide
+println("  Polynomial:         ", round(pvalue(mc_poly); digits=3)) # hide
+println("  Gaussian (Custom):  ", round(pvalue(mc_gauss); digits=3)) # hide
+
+# A small p-value provides evidence against the null hypothesis that the model captures the
+# observed firing pattern, while a large p-value indicates the model is consistent with the data.
+
 # ## Model Comparison
 # Let's compare all models by computing their negative log-likelihoods:
 
diff --git a/src/univariate/poisson/inhomogeneous/inhomogeneous_poisson_process.jl b/src/univariate/poisson/inhomogeneous/inhomogeneous_poisson_process.jl
index 0ce8e44..c56992b 100644
--- a/src/univariate/poisson/inhomogeneous/inhomogeneous_poisson_process.jl
+++ b/src/univariate/poisson/inhomogeneous/inhomogeneous_poisson_process.jl
@@ -107,6 +107,27 @@ function integrated_ground_intensity(
     )
 end
 
+## Time change
+
+"""
+    time_change(h::History, pp::InhomogeneousPoissonProcess)
+
+Apply the time rescaling theorem to history `h` using the compensator
+``Λ(t) = ∫_{tmin}^{t} λ(s) ds`` of the inhomogeneous Poisson process `pp`.
+
+The transformed event times form a unit-rate (homogeneous) Poisson process on
+``[0, Λ(tmax)]``, which makes the standard goodness-of-fit tests (e.g. KS against
+`Uniform` or `Exponential`) applicable.
+"""
+function time_change(h::History, pp::InhomogeneousPoissonProcess)
+    f = pp.intensity_function
+    config = pp.integration_config
+    new_tmax = integrated_intensity(f, h.tmin, h.tmax, config)
+    T = typeof(new_tmax)
+    new_times = T[integrated_intensity(f, h.tmin, t, config) for t in h.times]
+    return History(; times=new_times, tmin=zero(T), tmax=new_tmax, marks=h.marks)
+end
+
 ## Simulation
 function simulate(rng::AbstractRNG, pp::InhomogeneousPoissonProcess, tmin::Real, tmax::Real)
     return simulate_ogata(rng, pp, tmin, tmax)
diff --git a/test/hypothesis_tests.jl b/test/hypothesis_tests.jl
index 76cfa9e..e64eecd 100644
--- a/test/hypothesis_tests.jl
+++ b/test/hypothesis_tests.jl
@@ -68,3 +68,40 @@ end
 
     @test all(sort(test1.sim_stats) .== sort(test2.sim_stats)) # Test reproducibility
 end
+
+@testset "Inhomogeneous Poisson goodness-of-fit" begin
+    # KSDistance with InhomogeneousPoissonProcess relies on time_change.
+    # These tests pin down that the integration-based time_change feeds correctly
+    # through statistic / MonteCarloTest.
+
+    rng = Random.seed!(31415)
+    intensity_true = PolynomialIntensity([2.0, 0.3])
+    pp_true = InhomogeneousPoissonProcess(intensity_true)
+    h = simulate(rng, pp_true, 0.0, 50.0)
+
+    @testset "Statistic computation" begin
+        s_unif = statistic(KSDistance{Uniform}, pp_true, h)
+        s_exp = statistic(KSDistance{Exponential}, pp_true, h)
+        @test 0.0 <= s_unif <= 1.0
+        @test 0.0 <= s_exp <= 1.0
+    end
+
+    @testset "MonteCarloTest with correct model" begin
+        # Under the true model, the KS test should not reject at conventional levels.
+        mc = MonteCarloTest(
+            KSDistance{Exponential}, pp_true, h; n_sims=200, rng=Random.seed!(7)
+        )
+        @test mc.n_sims == 200
+        @test pvalue(mc) > 0.05
+    end
+
+    @testset "MonteCarloTest with misspecified model" begin
+        # A wildly wrong rate should make the time-rescaled times poorly Exp(1)
+        # under that model, so the test should reject.
+        pp_wrong = InhomogeneousPoissonProcess(PolynomialIntensity([20.0]))  # constant 20.0
+        mc = MonteCarloTest(
+            KSDistance{Exponential}, pp_wrong, h; n_sims=200, rng=Random.seed!(7)
+        )
+        @test pvalue(mc) < 0.05
+    end
+end
diff --git a/test/inhomogeneous_poisson_process.jl b/test/inhomogeneous_poisson_process.jl
index 9d18257..fea3f7b 100644
--- a/test/inhomogeneous_poisson_process.jl
+++ b/test/inhomogeneous_poisson_process.jl
@@ -9,6 +9,76 @@ using Test
 
 rng = Random.seed!(12345)
 
+@testset "Time change" begin
+    @testset "Polynomial intensity" begin
+        # λ(t) = 1 + 0.5*t  =>  Λ(t) - Λ(0) = t + 0.25*t²
+        intensity_linear = PolynomialIntensity([1.0, 0.5])
+        pp = InhomogeneousPoissonProcess(intensity_linear, Normal())
+        h = History([1.0, 2.0, 4.0], 0.0, 5.0, [0.1, 0.2, 0.3])
+
+        h_transf = time_change(h, pp)
+        expected_times = [t + 0.25 * t^2 for t in h.times]
+        expected_tmax = 5.0 + 0.25 * 25.0
+
+        @test h_transf.tmin ≈ 0.0
+        @test h_transf.tmax ≈ expected_tmax rtol = 1e-6
+        @test h_transf.times ≈ expected_times rtol = 1e-6
+        @test h_transf.marks == h.marks
+    end
+
+    @testset "Polynomial intensity with non-zero tmin" begin
+        # λ(t) = 2.0  =>  Λ(t) - Λ(tmin) = 2*(t - tmin)
+        intensity_const = PolynomialIntensity([2.0])
+        pp = InhomogeneousPoissonProcess(intensity_const)
+        h = History([3.0, 5.0, 7.0], 2.0, 10.0)
+
+        h_transf = time_change(h, pp)
+        @test h_transf.tmin ≈ 0.0
+        @test h_transf.tmax ≈ 2.0 * (10.0 - 2.0) rtol = 1e-6
+        @test h_transf.times ≈ [2.0 * (t - 2.0) for t in h.times] rtol = 1e-6
+    end
+
+    @testset "Empty history" begin
+        intensity_linear = PolynomialIntensity([1.0, 0.5])
+        pp = InhomogeneousPoissonProcess(intensity_linear)
+        h_empty = History(Float64[], 0.0, 5.0)
+
+        h_transf = time_change(h_empty, pp)
+        @test isempty(h_transf.times)
+        @test h_transf.tmin ≈ 0.0
+        @test h_transf.tmax ≈ 5.0 + 0.25 * 25.0 rtol = 1e-6
+    end
+
+    @testset "Custom intensity function" begin
+        # λ(t) = 1 + sin(t)^2   (no analytical integral, uses numerical)
+        custom_func = t -> 1.0 + sin(t)^2
+        pp = InhomogeneousPoissonProcess(custom_func)
+        h = History([0.5, 1.0, 2.0], 0.0, 3.0)
+
+        h_transf = time_change(h, pp)
+        @test h_transf.tmin ≈ 0.0
+        @test issorted(h_transf.times)
+        # Λ(t) for 1 + sin²(t) = t + (t - sin(t)cos(t))/2 = 1.5*t - 0.5*sin(t)*cos(t)
+        expected_times = [1.5 * t - 0.5 * sin(t) * cos(t) for t in h.times]
+        expected_tmax = 1.5 * 3.0 - 0.5 * sin(3.0) * cos(3.0)
+        @test h_transf.times ≈ expected_times rtol = 1e-6
+        @test h_transf.tmax ≈ expected_tmax rtol = 1e-6
+    end
+
+    @testset "Goodness-of-fit consistency" begin
+        # Simulating from a process and rescaling should yield approximately
+        # a unit-rate Poisson process: KS test against Uniform should not reject.
+        intensity_true = PolynomialIntensity([2.0, 0.3])
+        pp_true = InhomogeneousPoissonProcess(intensity_true, Dirac(nothing))
+        h = simulate(rng, pp_true, 0.0, 50.0)
+
+        h_transf = time_change(h, pp_true)
+        @test h_transf.tmin ≈ 0.0
+        @test issorted(h_transf.times)
+        @test all(h_transf.tmin .<= h_transf.times .<= h_transf.tmax)
+    end
+end
+
 @testset "PolynomialIntensity" begin
     # Linear intensity: λ(t) = 1 + 0.5*t
     intensity_linear = PolynomialIntensity([1.0, 0.5])

From 8635d077dcd8e9724218dc879357500c31ae6842 Mon Sep 17 00:00:00 2001
From: Ryan Senne <50930199+rsenne@users.noreply.github.com>
Date: Fri, 8 May 2026 23:09:27 -0400
Subject: [PATCH 2/6] Fixes to copilot issues

---
 .../inhomogeneous_poisson_process.jl          | 15 +++++++++++-
 test/hypothesis_tests.jl                      | 24 ++++++++++---------
 test/inhomogeneous_poisson_process.jl         |  2 --
 3 files changed, 27 insertions(+), 14 deletions(-)

diff --git a/src/univariate/poisson/inhomogeneous/inhomogeneous_poisson_process.jl b/src/univariate/poisson/inhomogeneous/inhomogeneous_poisson_process.jl
index c56992b..bf3962d 100644
--- a/src/univariate/poisson/inhomogeneous/inhomogeneous_poisson_process.jl
+++ b/src/univariate/poisson/inhomogeneous/inhomogeneous_poisson_process.jl
@@ -125,7 +125,20 @@ function time_change(h::History, pp::InhomogeneousPoissonProcess)
     new_tmax = integrated_intensity(f, h.tmin, h.tmax, config)
     T = typeof(new_tmax)
     new_times = T[integrated_intensity(f, h.tmin, t, config) for t in h.times]
-    return History(; times=new_times, tmin=zero(T), tmax=new_tmax, marks=h.marks)
+    # `new_tmax` and each `new_times[i]` come from independent quadrature calls,
+    # so solver / floating-point error can let a transformed time slightly exceed
+    # `new_tmax`. Widen `new_tmax` to the observed maximum and skip History's
+    # bounds check so events are never silently dropped.
+    if !isempty(new_times)
+        new_tmax = max(new_tmax, last(new_times) * (one(T) + eps(T)))
+    end
+    return History(;
+        times=new_times,
+        tmin=zero(T),
+        tmax=new_tmax,
+        marks=h.marks,
+        check_args=false,
+    )
 end
 
 ## Simulation
diff --git a/test/hypothesis_tests.jl b/test/hypothesis_tests.jl
index e64eecd..f640f91 100644
--- a/test/hypothesis_tests.jl
+++ b/test/hypothesis_tests.jl
@@ -86,22 +86,24 @@ end
         @test 0.0 <= s_exp <= 1.0
     end
 
-    @testset "MonteCarloTest with correct model" begin
-        # Under the true model, the KS test should not reject at conventional levels.
-        mc = MonteCarloTest(
+    @testset "MonteCarloTest correct vs. misspecified" begin
+        # correct model should fit the data substantially better than a wildly
+        # misspecified one.
+        mc_true = MonteCarloTest(
             KSDistance{Exponential}, pp_true, h; n_sims=200, rng=Random.seed!(7)
         )
-        @test mc.n_sims == 200
-        @test pvalue(mc) > 0.05
-    end
+        @test mc_true.n_sims == 200
 
-    @testset "MonteCarloTest with misspecified model" begin
-        # A wildly wrong rate should make the time-rescaled times poorly Exp(1)
-        # under that model, so the test should reject.
         pp_wrong = InhomogeneousPoissonProcess(PolynomialIntensity([20.0]))  # constant 20.0
-        mc = MonteCarloTest(
+        mc_wrong = MonteCarloTest(
             KSDistance{Exponential}, pp_wrong, h; n_sims=200, rng=Random.seed!(7)
         )
-        @test pvalue(mc) < 0.05
+
+        # The misspecified model's KS statistic should be substantially
+        # larger than the correct one's — this is the most direct measure
+        # of relative fit quality and is robust to specific p-value seeds.
+        @test mc_wrong.stat > 2 * mc_true.stat
+        # Correct model's p-value should also exceed the misspecified one's.
+        @test pvalue(mc_true) > pvalue(mc_wrong)
     end
 end
diff --git a/test/inhomogeneous_poisson_process.jl b/test/inhomogeneous_poisson_process.jl
index fea3f7b..efc3d12 100644
--- a/test/inhomogeneous_poisson_process.jl
+++ b/test/inhomogeneous_poisson_process.jl
@@ -50,7 +50,6 @@ rng = Random.seed!(12345)
     end
 
     @testset "Custom intensity function" begin
-        # λ(t) = 1 + sin(t)^2   (no analytical integral, uses numerical)
         custom_func = t -> 1.0 + sin(t)^2
         pp = InhomogeneousPoissonProcess(custom_func)
         h = History([0.5, 1.0, 2.0], 0.0, 3.0)
@@ -58,7 +57,6 @@ rng = Random.seed!(12345)
         h_transf = time_change(h, pp)
         @test h_transf.tmin ≈ 0.0
         @test issorted(h_transf.times)
-        # Λ(t) for 1 + sin²(t) = t + (t - sin(t)cos(t))/2 = 1.5*t - 0.5*sin(t)*cos(t)
         expected_times = [1.5 * t - 0.5 * sin(t) * cos(t) for t in h.times]
         expected_tmax = 1.5 * 3.0 - 0.5 * sin(3.0) * cos(3.0)
         @test h_transf.times ≈ expected_times rtol = 1e-6

From 6d7814eb2258afc440f1e13535b5e7cae14a0450 Mon Sep 17 00:00:00 2001
From: Ryan Senne <50930199+rsenne@users.noreply.github.com>
Date: Fri, 8 May 2026 23:27:50 -0400
Subject: [PATCH 3/6] Update
 src/univariate/poisson/inhomogeneous/inhomogeneous_poisson_process.jl

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>
---
 .../poisson/inhomogeneous/inhomogeneous_poisson_process.jl  | 6 +-----
 1 file changed, 1 insertion(+), 5 deletions(-)

diff --git a/src/univariate/poisson/inhomogeneous/inhomogeneous_poisson_process.jl b/src/univariate/poisson/inhomogeneous/inhomogeneous_poisson_process.jl
index bf3962d..261e2c8 100644
--- a/src/univariate/poisson/inhomogeneous/inhomogeneous_poisson_process.jl
+++ b/src/univariate/poisson/inhomogeneous/inhomogeneous_poisson_process.jl
@@ -133,11 +133,7 @@ function time_change(h::History, pp::InhomogeneousPoissonProcess)
         new_tmax = max(new_tmax, last(new_times) * (one(T) + eps(T)))
     end
     return History(;
-        times=new_times,
-        tmin=zero(T),
-        tmax=new_tmax,
-        marks=h.marks,
-        check_args=false,
+        times=new_times, tmin=zero(T), tmax=new_tmax, marks=h.marks, check_args=false
     )
 end
 

From 4bb42d8bf30ba85f3a18ad6e05b83ff858b5d22b Mon Sep 17 00:00:00 2001
From: Ryan Senne <50930199+rsenne@users.noreply.github.com>
Date: Sat, 9 May 2026 09:27:21 -0400
Subject: [PATCH 4/6] Workaround for bug in histroy

---
 .../inhomogeneous/inhomogeneous_poisson_process.jl | 14 +++++++-------
 test/runtests.jl                                   |  4 ++--
 2 files changed, 9 insertions(+), 9 deletions(-)

diff --git a/src/univariate/poisson/inhomogeneous/inhomogeneous_poisson_process.jl b/src/univariate/poisson/inhomogeneous/inhomogeneous_poisson_process.jl
index 261e2c8..bf87dab 100644
--- a/src/univariate/poisson/inhomogeneous/inhomogeneous_poisson_process.jl
+++ b/src/univariate/poisson/inhomogeneous/inhomogeneous_poisson_process.jl
@@ -126,15 +126,15 @@ function time_change(h::History, pp::InhomogeneousPoissonProcess)
     T = typeof(new_tmax)
     new_times = T[integrated_intensity(f, h.tmin, t, config) for t in h.times]
     # `new_tmax` and each `new_times[i]` come from independent quadrature calls,
-    # so solver / floating-point error can let a transformed time slightly exceed
-    # `new_tmax`. Widen `new_tmax` to the observed maximum and skip History's
-    # bounds check so events are never silently dropped.
+    # so solver / floating-point error can both swap adjacent transformed
+    # events and let one of them slightly exceed `new_tmax`. Widen `new_tmax`
+    # to the observed maximum (with an `eps` cushion) so no event is dropped
+    # at the boundary; small order drift is then handled by the History
+    # constructor's automatic `sortperm` pass under `check_args=true`.
     if !isempty(new_times)
-        new_tmax = max(new_tmax, last(new_times) * (one(T) + eps(T)))
+        new_tmax = max(new_tmax, maximum(new_times) * (one(T) + eps(T)))
     end
-    return History(;
-        times=new_times, tmin=zero(T), tmax=new_tmax, marks=h.marks, check_args=false
-    )
+    return History(; times=new_times, tmin=zero(T), tmax=new_tmax, marks=h.marks)
 end
 
 ## Simulation
diff --git a/test/runtests.jl b/test/runtests.jl
index 762da97..1203e04 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -20,8 +20,8 @@ DocMeta.setdocmeta!(PointProcesses, :DocTestSetup, :(using PointProcesses); recu
         Aqua.test_all(PointProcesses; ambiguities=false, deps_compat=(; check_extras=false))
     end
     @testset verbose = false "Code Linting" begin
-        # test on 1.11 at a minimum and not on pre-release 
-        if VERSION >= v"1.11" && isempty(VERSION.prerelease)
+        # Skip JET on Julia pre-releases (where JET often hasn't caught up yet).
+        if isempty(VERSION.prerelease)
             JET.test_package(PointProcesses; target_modules=(PointProcesses,))
         end
     end

From db3aadf9f0fe98d513616c56a8f26193be471f15 Mon Sep 17 00:00:00 2001
From: Ryan Senne <50930199+rsenne@users.noreply.github.com>
Date: Sat, 9 May 2026 10:09:05 -0400
Subject: [PATCH 5/6] Attempt at JET fix

---
 test/Project.toml |  1 -
 test/runtests.jl  | 11 +++++++++--
 2 files changed, 9 insertions(+), 3 deletions(-)

diff --git a/test/Project.toml b/test/Project.toml
index a8f20a9..23ebec7 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -4,7 +4,6 @@ DensityInterface = "b429d917-457f-4dbc-8f4c-0cc954292b1d"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
-JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
diff --git a/test/runtests.jl b/test/runtests.jl
index 1203e04..c5e41ec 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -3,8 +3,8 @@ using DensityInterface
 using Documenter
 using Distributions
 using ForwardDiff
-using JET
 using LinearAlgebra
+using Pkg
 using PointProcesses
 using Random
 using Statistics
@@ -20,8 +20,15 @@ DocMeta.setdocmeta!(PointProcesses, :DocTestSetup, :(using PointProcesses); recu
         Aqua.test_all(PointProcesses; ambiguities=false, deps_compat=(; check_extras=false))
     end
     @testset verbose = false "Code Linting" begin
-        # Skip JET on Julia pre-releases (where JET often hasn't caught up yet).
+        # Skip JET on Julia pre-releases (where JET typically hasn't caught up
+        # yet and there is no compatible JET version to resolve against). JET is
+        # not listed in test/Project.toml's [deps] for the same reaso, having
+        # it there would make `Pkg.test` fail at the resolution step on
+        # prereleases, before this guard ever runs. Install on demand only when
+        # we're on a stable Julia.
         if isempty(VERSION.prerelease)
+            Pkg.add("JET")
+            @eval using JET
             JET.test_package(PointProcesses; target_modules=(PointProcesses,))
         end
     end

From 77740645c059102a745d7ea80c8fd3857bb2d667 Mon Sep 17 00:00:00 2001
From: Ryan Senne <50930199+rsenne@users.noreply.github.com>
Date: Sat, 9 May 2026 10:11:55 -0400
Subject: [PATCH 6/6] Add Pkg to requiremtns

---
 test/Project.toml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/test/Project.toml b/test/Project.toml
index 23ebec7..98b3491 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -6,6 +6,7 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
+Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsAPI = "82ae8749-77ed-4fe6-ae5f-f523153014b0"