Julia implementations for `poiscdf` and `poislogcdf` #112

torfjelde · 2021-04-06T18:13:06Z

This PR adds Julia implementations for poiscdf and poislogcdf.

I'm a bit uncertain where tests for these should go, so I would appreciate some advice on that:)

codecov-io · 2021-04-06T18:15:14Z

Codecov Report

Merging #112 (a4a58dc) into master (bc45e18) will decrease coverage by 0.24%.
The diff coverage is 0.00%.

@@            Coverage Diff             @@
##           master     #112      +/-   ##
==========================================
- Coverage   39.74%   39.50%   -0.25%     
==========================================
  Files          12       12              
  Lines         478      481       +3     
==========================================
  Hits          190      190              
- Misses        288      291       +3

Impacted Files	Coverage Δ
src/distrs/pois.jl	`50.00% <0.00%> (-50.00%)`	⬇️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update bc45e18...a4a58dc. Read the comment docs.

src/distrs/pois.jl

devmotion · 2021-04-06T19:14:55Z

The generic implementations of (log)pdf are compared with Rmath in https://github.com/JuliaStats/StatsFuns.jl/blob/master/test/generic.jl. Similarly, one should check the implementations of (log)cdf.

Co-authored-by: David Widmann <devmotion@users.noreply.github.com>

torfjelde · 2021-04-06T19:23:46Z

The generic implementations of (log)pdf are compared with Rmath in https://github.com/JuliaStats/StatsFuns.jl/blob/master/test/generic.jl. Similarly, one should check the implementations of (log)cdf.

Cool 👍 But there are some that already have these implemented, e.g. normcdf, but there doesn't exist a test similar to genericcomp for cdf already?

EDIT: If that's not the case, I'll imlement one ofc. I was just making sure given that there are already similar methods which I assumed were already tested.

devmotion · 2021-04-06T19:24:14Z

No, I think it is missing.

Co-authored-by: David Widmann <devmotion@users.noreply.github.com>

torfjelde · 2021-04-07T12:31:21Z

src/distrs/pois.jl

+poiscdf(λ::Number, x::Number; precision=0) = last(gamma_inc(floor(x + 1), λ, precision))
+poisccdf(λ::Number, x::Number; precision=0) = first(gamma_inc(floor(x + 1), λ, precision))


Btw, should precision be argument rather than kwarg? @devmotion

It doesn't matter but IMO it is a bit more natural as a keyword argument since we don't dispatch on it.

I don't think we should expose the precision argument here. The implementation in SpecialFunctions happens to support it but it's not something that is generally available I don't think it will be in demand here. If it is then I think it will be more natural to get the lower precision version as part of a method for Float32.

torfjelde · 2021-04-07T12:37:34Z

src/distrs/pois.jl

@@ -19,6 +19,13 @@ poislogpdf(λ::T, x::T) where {T <: Real} = xlogy(x, λ) - λ - loggamma(x + 1)

 poislogpdf(λ::Number, x::Number) = poislogpdf(promote(float(λ), x)...)


I think I get why this is done btw: loggamma(::ForwardDiff.Dual, ::Int) won't have a diff-rule defined for it.

Actually, it doesn't seem like these have rules defined in DiffRules.jl!

There is no DiffRule defined for loggamma and gamma_inc anyway so it shouldn't matter at the moment 🤷 I added ChainRules definitions to SpecialFunctions (JuliaMath/SpecialFunctions.jl#305) but currently they can't be used with Zygote due to FluxML/Zygote.jl#873.

So can we avoid these conversions if possible for now and only add them later if needed?

I'll have to look at it later, but it seems like there are particular execution branches (depending on the values) which ForwardDiff.jl fails to automatically obtain derivatives for. So when I just tried a couple of values everything seemed to work, but when I let HMC loose on a distribution using this function under the hood, it eventually runs into a MethodError.

Haha, man, you're always like one step ahead: I just looked at your PR and that's exactly the "solution" i arrived at literally 10mins ago, i.e. we can only take derivative wrt. second argument but that's also what we care about so w/e.

There is no DiffRule defined for loggamma and gamma_inc anyway so it shouldn't matter at the moment

But wouldn't adding those fix this issue? Same as you've done in ChainRules now?

Sure, if I used the correct derivatives it should be fine 😛 AFAICT the convention in DiffRules would be to set the derivatives of the first argument to NaN, if it is not implemented properly.

andreasnoack · 2021-05-13T20:01:34Z

Regarding the testing then see #113. I think we should just start to use the native implementations unconditionally. Then they'll be tested automatically against the Rmath versions with the current test setup.

andreasnoack · 2021-05-13T20:10:14Z

src/distrs/pois.jl

@@ -19,6 +19,16 @@ poislogpdf(λ::T, x::T) where {T <: Real} = xlogy(x, λ) - λ - loggamma(x + 1)

 poislogpdf(λ::Number, x::Number) = poislogpdf(promote(float(λ), x)...)

+function poislogcdf(λ::T, x::T) where {T<:Real}
+    return loggamma(floor(x + 1), float(λ)) - logfactorial(floor(x))


Very minor detail, but it seems a little odd to me use logfactorial here together with loggamma. In particular when we also use loggamma(x+1) in the pdf.

andreasnoack · 2021-05-17T19:47:02Z

@torfjelde It would be great to have this one completed. Do you expect to have the cycles to complete it soon?

torfjelde · 2021-05-21T09:34:34Z

It would be great to have this one completed. Do you expect to have the cycles to complete it soon?

Sorry, been quite busy lately. I'm travelling this weekend so won't be able to do much, but I'll try to get it done early next week if that's okay?

andreasnoack · 2021-05-21T10:20:48Z

No worries. I ended up doing it as part of #113

torfjelde · 2021-06-07T07:41:43Z

Btw, should we just close this PR then?

andreasnoack · 2021-06-07T08:00:04Z

Yeah. Let's do that (and then I'll have to get #113 completed)

torfjelde added 2 commits April 6, 2021 20:10

added Julia implementations of poislogcdf and poiscdf

199ea18

ensure λ is float in poislogcdf

1938e52

devmotion reviewed Apr 6, 2021

View reviewed changes

src/distrs/pois.jl Outdated Show resolved Hide resolved

devmotion reviewed Apr 6, 2021

View reviewed changes

src/distrs/pois.jl Outdated Show resolved Hide resolved

devmotion reviewed Apr 6, 2021

View reviewed changes

src/distrs/pois.jl Outdated Show resolved Hide resolved

Apply suggestions from code review

a4a58dc

Co-authored-by: David Widmann <devmotion@users.noreply.github.com>

Update src/distrs/pois.jl

ae4391e

Co-authored-by: David Widmann <devmotion@users.noreply.github.com>

torfjelde commented Apr 7, 2021

View reviewed changes

should not dispatch correctly for diffrules

0bae7bf

andreasnoack mentioned this pull request May 13, 2021

Use julia implementations for pdfs and some cdf-like functions #113

Merged

andreasnoack reviewed May 13, 2021

View reviewed changes

andreasnoack closed this Jun 7, 2021

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Julia implementations for `poiscdf` and `poislogcdf` #112

Julia implementations for `poiscdf` and `poislogcdf` #112

torfjelde commented Apr 6, 2021

codecov-io commented Apr 6, 2021 •

edited

devmotion commented Apr 6, 2021

torfjelde commented Apr 6, 2021 •

edited

devmotion commented Apr 6, 2021

torfjelde Apr 7, 2021

devmotion Apr 7, 2021

andreasnoack May 13, 2021

torfjelde Apr 7, 2021

torfjelde Apr 7, 2021

devmotion Apr 7, 2021

torfjelde Apr 7, 2021

torfjelde Apr 7, 2021

devmotion Apr 7, 2021

andreasnoack commented May 13, 2021

andreasnoack May 13, 2021

andreasnoack commented May 17, 2021

torfjelde commented May 21, 2021

andreasnoack commented May 21, 2021

torfjelde commented Jun 7, 2021 •

edited

andreasnoack commented Jun 7, 2021

		poiscdf(λ::Number, x::Number; precision=0) = last(gamma_inc(floor(x + 1), λ, precision))
		poisccdf(λ::Number, x::Number; precision=0) = first(gamma_inc(floor(x + 1), λ, precision))

		@@ -19,6 +19,13 @@ poislogpdf(λ::T, x::T) where {T <: Real} = xlogy(x, λ) - λ - loggamma(x + 1)

		poislogpdf(λ::Number, x::Number) = poislogpdf(promote(float(λ), x)...)

Julia implementations for poiscdf and poislogcdf #112

Julia implementations for poiscdf and poislogcdf #112

Conversation

torfjelde commented Apr 6, 2021

codecov-io commented Apr 6, 2021 • edited

Codecov Report

devmotion commented Apr 6, 2021

torfjelde commented Apr 6, 2021 • edited

devmotion commented Apr 6, 2021

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

andreasnoack commented May 13, 2021

Choose a reason for hiding this comment

andreasnoack commented May 17, 2021

torfjelde commented May 21, 2021

andreasnoack commented May 21, 2021

torfjelde commented Jun 7, 2021 • edited

andreasnoack commented Jun 7, 2021

Julia implementations for `poiscdf` and `poislogcdf` #112

Julia implementations for `poiscdf` and `poislogcdf` #112

codecov-io commented Apr 6, 2021 •

edited

torfjelde commented Apr 6, 2021 •

edited

torfjelde commented Jun 7, 2021 •

edited