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..bf87dab 100644 --- a/src/univariate/poisson/inhomogeneous/inhomogeneous_poisson_process.jl +++ b/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) diff --git a/test/Project.toml b/test/Project.toml index a8f20a9..98b3491 100644 --- a/test/Project.toml +++ b/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" diff --git a/test/hypothesis_tests.jl b/test/hypothesis_tests.jl index 76cfa9e..f640f91 100644 --- a/test/hypothesis_tests.jl +++ b/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 diff --git a/test/inhomogeneous_poisson_process.jl b/test/inhomogeneous_poisson_process.jl index 9d18257..efc3d12 100644 --- a/test/inhomogeneous_poisson_process.jl +++ b/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]) diff --git a/test/runtests.jl b/test/runtests.jl index 762da97..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 - # 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