Add rand safe version of Poisson and Negative binomial distributions
#418
Conversation
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Try this Pull Request!Open Julia and type: import Pkg
Pkg.activate(temp=true)
Pkg.add(url="https://github.com/CDCgov/Rt-without-renewal", rev="safe-discrete-dists", subdir="EpiAware")
using EpiAware |
Codecov ReportAttention: Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## remove-nan-handling #418 +/- ##
=======================================================
- Coverage 93.71% 88.67% -5.05%
=======================================================
Files 54 56 +2
Lines 557 715 +158
=======================================================
+ Hits 522 634 +112
- Misses 35 81 +46 ☔ View full report in Codecov by Sentry. |
Benchmark resultJudge resultBenchmark Report for /home/runner/work/Rt-without-renewal/Rt-without-renewalJob Properties
ResultsA ratio greater than
Benchmark Group ListHere's a list of all the benchmark groups executed by this job:
Julia versioninfoTargetBaselineTarget resultBenchmark Report for /home/runner/work/Rt-without-renewal/Rt-without-renewalJob Properties
ResultsBelow is a table of this job's results, obtained by running the benchmarks.
Benchmark Group ListHere's a list of all the benchmark groups executed by this job:
Julia versioninfoBaseline resultBenchmark Report for /home/runner/work/Rt-without-renewal/Rt-without-renewalJob Properties
ResultsBelow is a table of this job's results, obtained by running the benchmarks.
Benchmark Group ListHere's a list of all the benchmark groups executed by this job:
Julia versioninfoRuntime information
|
|
Note that this is branched from We can either address that with a new issue, or use that as an opportunity to test using the new dists? |
seabbs
reviewed
Aug 14, 2024
There was a problem hiding this comment.
I think this looks good and I think we should replace our usage everywhere with this (so that i.e Documenter passes). I also think we should push ideas about reducing copy and paste here into a issue but definitely think it should be investigated as it should help with our partial pooling problems
seabbs
reviewed
Aug 14, 2024
|
This needs a benchmark before we merge it. |
If we call |
|
I got irritated by a stochastic CI failure on how different the sample variance is from sampling from |
Benchmark resultJudge resultBenchmark Report for /home/runner/work/Rt-without-renewal/Rt-without-renewalJob Properties
ResultsA ratio greater than
Benchmark Group ListHere's a list of all the benchmark groups executed by this job:
Julia versioninfoTargetBaselineTarget resultBenchmark Report for /home/runner/work/Rt-without-renewal/Rt-without-renewalJob Properties
ResultsBelow is a table of this job's results, obtained by running the benchmarks.
Benchmark Group ListHere's a list of all the benchmark groups executed by this job:
Julia versioninfoBaseline resultBenchmark Report for /home/runner/work/Rt-without-renewal/Rt-without-renewalJob Properties
ResultsBelow is a table of this job's results, obtained by running the benchmarks.
Benchmark Group ListHere's a list of all the benchmark groups executed by this job:
Julia versioninfoRuntime information
|
Benchmark resultJudge resultBenchmark Report for /home/runner/work/Rt-without-renewal/Rt-without-renewalJob Properties
ResultsA ratio greater than
Benchmark Group ListHere's a list of all the benchmark groups executed by this job:
Julia versioninfoTargetBaselineTarget resultBenchmark Report for /home/runner/work/Rt-without-renewal/Rt-without-renewalJob Properties
ResultsBelow is a table of this job's results, obtained by running the benchmarks.
Benchmark Group ListHere's a list of all the benchmark groups executed by this job:
Julia versioninfoBaseline resultBenchmark Report for /home/runner/work/Rt-without-renewal/Rt-without-renewalJob Properties
ResultsBelow is a table of this job's results, obtained by running the benchmarks.
Benchmark Group ListHere's a list of all the benchmark groups executed by this job:
Julia versioninfoRuntime information
|
|
Documenter CI works now. To recap:
|
Benchmark resultJudge resultBenchmark Report for /home/runner/work/Rt-without-renewal/Rt-without-renewalJob Properties
ResultsA ratio greater than
Benchmark Group ListHere's a list of all the benchmark groups executed by this job:
Julia versioninfoTargetBaselineTarget resultBenchmark Report for /home/runner/work/Rt-without-renewal/Rt-without-renewalJob Properties
ResultsBelow is a table of this job's results, obtained by running the benchmarks.
Benchmark Group ListHere's a list of all the benchmark groups executed by this job:
Julia versioninfoBaseline resultBenchmark Report for /home/runner/work/Rt-without-renewal/Rt-without-renewalJob Properties
ResultsBelow is a table of this job's results, obtained by running the benchmarks.
Benchmark Group ListHere's a list of all the benchmark groups executed by this job:
Julia versioninfoRuntime information
|
seabbs
reviewed
Aug 15, 2024
| p = μ / (μ + ex_σ²) | ||
| r = μ² / ex_σ² | ||
| return NegativeBinomial(r, p) | ||
| σ² = μ + α * μ² |
seabbs
approved these changes
Aug 15, 2024
This was referenced Aug 15, 2024
github-merge-queue bot
pushed a commit
that referenced
this pull request
Aug 15, 2024
* remove nan handling * get rid of all NB padding and clamping * remove nan handling * get rid of all NB padding and clamping * remove overflow test * Add `rand` safe version of Poisson and Negative binomial distributions (#418) * SafePoisson with safety for large means * better selection for conversion to Int or BigInt * add SafeNegativeBinomial * add unit tests to doctests * reformat * Add type promotion so AD works with distribution constructor * Add logpdf grad call unit tests for Safe discrete dists * reformat * change neg bin param to (r, p) * Update utils.jl * reformat * change empirical var test to more principled approach * add default nadapts rather than just 50% of target sampling * Update NUTSampler.jl * set dist check_args = false * Set nadapts to Turing Default * reformat --------- Co-authored-by: Samuel Brand <[email protected]>
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
This PR contributes a possible solution to #415 by creating distributions that work with
Distributionsbut don't trigger aInexactErrorwhenrandis called on them with a large mean value (noting that a large mean value can triggered by random walks on the log-mean). These distributions are:SafePoissonSafeNegativeBinomialLargely this was boilerplate (if we do more of this we should consider a macro to avoid this), with the difference that I replace
floor(Int, x)calls with a helper function that determines whetherfloor(BigInt, x)should be used rather than throwing anInexactError.Notable features/limitations
PoissonRandomafter checking that its still faster thanDistribution.Poissonsampling.SafeNegativeBinomialI'm using the parameterisation with meanMy reasoning here is in the doc string; I think this is the easiest parameterisation to reason on priors because$\alpha \approx \text{std} / \text{mean}$ when the mean is larger than $1/ \alpha$ . This differs from
NegativeBinomialMeanClustso I'm happy to change.ReverseDiffandForwardDiffto the test env. This was because we need to be sure that we can backprop throughlogpdfcalls withSafePoissonorSafeNegativeBinomialdistribution structs. I've added unit tests here with very large mean values to check validity in extreme parameter regimes.randcalls are slower to these distributions (but not done large formal analysis).NegativeBinomialMeanClust.I know @damonbayer has been interested in Neg bin safety too.