Add time rescaling to inhomogeneous

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:
 

src/univariate/poisson/inhomogeneous/inhomogeneous_poisson_process.jl

@@ -107,6 +107,36 @@ 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]
+    # `new_tmax` and each `new_times[i]` come from independent quadrature calls,
+    # 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, maximum(new_times) * (one(T) + eps(T)))
+    end
+    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)

test/Project.toml

@@ -4,9 +4,9 @@ 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"
+Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsAPI = "82ae8749-77ed-4fe6-ae5f-f523153014b0"

test/hypothesis_tests.jl

@@ -68,3 +68,42 @@ 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 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_true.n_sims == 200
+
+        pp_wrong = InhomogeneousPoissonProcess(PolynomialIntensity([20.0]))  # constant 20.0
+        mc_wrong = MonteCarloTest(
+            KSDistance{Exponential}, pp_wrong, h; n_sims=200, rng=Random.seed!(7)
+        )
+
+        # 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

test/inhomogeneous_poisson_process.jl

@@ -9,6 +9,74 @@ 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
+        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)
+        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])

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
-        # 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 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