Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add lbinomial function based on lbeta #99

Merged
merged 2 commits into from
Jul 23, 2018

Conversation

freeboson
Copy link
Contributor

This was originally submitted to base in JuliaLang/julia#27419, but JuliaLang/julia#27473 moved special/gamma.jl things here, so I'm making the same PR here now. Details are copied over:

This implementation of lbinomial uses the fact that binomial(n,k) = 1/( (n+1) * B(n - k + 1, k + 1) ). I have added tests for some simple cases, and then also comparing for very large n with binomial(::BigInt, ::BigInt).

The main motivation is evaluating log(binomial(n,k)) for large n with k near n/2 without approximating with Poisson:

julia> log(binomial(70, 25)) # OK
43.311504703669215

julia> log(binomial(70, 26)) # :(
ERROR: InexactError()
Stacktrace:
 [1] binomial(::Int64, ::Int64) at ./intfuncs.jl:810

julia> log(binomial(BigInt(70), BigInt(26))) # OK, if you can afford it
4.386007065541805521142055198415251198714576775237582336220765505130378175238536e+01

julia> ans - lbinomial(70, 26) # Close enough?
-8.056495743460440321531836165841374176637792344948696218247614638531207882439305e-15

In R, this function is available as lchoose(n, k) and is a builtin.

src/gamma.jl Outdated
S = float(T)
(k < 0) && return typemin(S)
if n < 0
n = -n + k - 1
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Incorrect indentation here (needs one more space)

src/gamma.jl Outdated
(k == 0 || k == n) && return zero(S)
(k == 1) && return log(abs(n))
if k > (n>>1)
k = n - k
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here too

src/gamma.jl Outdated
end
lbinomial(n::Integer, k::Integer) = lbinomial(promote(n, k)...)

export lbinomial
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should go with the other exports.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I put it now with cosint. Not sure if this is the "normal" place or not.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, that's fine.

test/runtests.jl Outdated
@@ -636,3 +636,17 @@ end

@test beta(big(1.0),big(1.2)) ≈ beta(1.0,1.2) rtol=4*eps()
end

@testset "lbinomial" begin
lbfixed(n) = k -> lbinomial(n, k)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why is this necessary? You can write lbfixed(200).(0:200) below as lbinomial.(200, 0:200).

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Note that the same applies to the function defined directly below this line.

@freeboson
Copy link
Contributor Author

Rebased with the corrections, dropped the lambdas in testing.

@simonbyrne simonbyrne merged commit a57c72e into JuliaMath:master Jul 23, 2018
@simonbyrne
Copy link
Member

Thanks!

@cossio
Copy link
Contributor

cossio commented Nov 4, 2023

Is it intentional that this function is not documented? Should it be considered public API?

@ararslan
Copy link
Member

ararslan commented Dec 7, 2023

It's deprecated in favor of logabsbinomial:

@deprecate lbinomial(x) logabsbinomial(x)[1]

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants