Add time rescaling to inhomogeneous#81
Conversation
Codecov Report✅ All modified and coverable lines are covered by tests. 📢 Thoughts on this report? Let us know! |
There was a problem hiding this comment.
Pull request overview
Adds a time-rescaling (“time change”) transform for InhomogeneousPoissonProcess so existing goodness-of-fit tests (KS vs Uniform / Exponential) can be applied after compensator-based rescaling, with accompanying tests and documentation updates.
Changes:
- Implement
time_change(h, pp::InhomogeneousPoissonProcess)using the integrated intensity (compensator). - Add unit tests for time rescaling on polynomial and custom intensities, plus integration into hypothesis-test workflows.
- Extend the inhomogeneous Poisson documentation example with a Monte Carlo KS goodness-of-fit section.
Reviewed changes
Copilot reviewed 4 out of 4 changed files in this pull request and generated 3 comments.
| File | Description |
|---|---|
src/univariate/poisson/inhomogeneous/inhomogeneous_poisson_process.jl |
Adds time_change for inhomogeneous Poisson processes via integrated_intensity. |
test/inhomogeneous_poisson_process.jl |
Adds direct tests for correctness of the time-change transform (incl. non-zero tmin, empty history, custom intensity). |
test/hypothesis_tests.jl |
Adds tests exercising KSDistance/MonteCarloTest integration for inhomogeneous Poisson via time_change. |
docs/examples/Inhomogeneous.jl |
Documents goodness-of-fit testing via time rescaling and Monte Carlo KS tests. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
…cess.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>
| 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 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 | ||
|
|
There was a problem hiding this comment.
If we have to check if the order was changed for new_tmax, then isn't it possible that the order of any pair of events got switched? In this case, should the History be initialized with check_args set to true?
There was a problem hiding this comment.
You are right. My mistake. After reviewing the code of History i think this speaks to an underlying issue with our logic for input handling. I think History might be better of as a validator and not a fixer. What i mean by this is I feel history might be better of only qualifying if the history times meet the proper checks (i.e., bounded properly and sorted).
My thinking for this is that there are three reasonable cases unsorted data happens: i.) a user has unsorted data and needs to sort it ii.) a user has an actual error in their data they don't know about iii) we have numerical issues that are introduced from independent quadrature calls.
Right now the behavior of History basically conflates all three as equally valid/equally reasonable occurrences. I don't think we should try to reason about what the case of the above is except when we know for sure what is happening like in numerical issue introduced by quadrature.
On our side my thought is we handle these numerical issues and then use History as a way to validate our own code works as expected (e.g., by sorting events and also expanding the tmax). In the case of the user i think we should just throw an error. That way they have to opt-in to any changes of their data and or reciprocally opt-out.
This is getting a bit off-topic for this specific case, but I think this is something that might be at the least thinking about more concretely. We could easily supply some easy helpers so that a user knows they are opting into these sort of data fixes e.g., sort_history! and restrict. This would be a breaking change but we aren't at major version 1 yet so semantic versioning doesn't require us to notify users but we could at least add a deprecation cycle if we wanted to.
Hope this makes sense--let me know if it doesn't! I can also move this over to a new issue.
There was a problem hiding this comment.
My answer to the History stuff will be in the main thread. But about the ordering, since we will need to call quadrature n+1 times anyway, can't we accumulate the integral of interarrival times instead of calculating the integral of each event time? I would guess quadrature would not return negative results for non negative functions, therefore this would not change the ordering.
|
Great. I was planning on implementing this next. |
|
I re-requested review as I made further changes. I edited the behavior per your comment (see the response I left i think it's worth thinking about), and also made JET runtime conditional as the constant failure on pre-release was annoying me lol. This PR got a little noisy as a result so I'd probably squash merge it |
|
About |
| # 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 |
|
On the cumulative inter-arrival idea: i think this is exactly right and we should do this in another PR. If we integrate over On keeping History as a fixer with a warning, i see the appeal, especially the worry about MC overhead, but im not certain this would cost us much. The actual validation is just Nonetheless, I'm going to merge this PR and I can start a discussion on these points as i view them as separate issues whether added or not |
Resolves #78