Add doc strings and explicit BigInt method for lstirling_asym. #36
Conversation
src/basicfuns.jl
Outdated
|
|
||
| Return `x * log(x)` for `x ≥ 0`. (Special case is `x = 0`.) |
There was a problem hiding this comment.
Maybe
Return `x * log(x)` for `x ≥ 0`, handling `x = 0` as the limit
```jldoctest
julia> StatsFuns.xlogx(0)
0.0
```
|
I've made a few minor tweaks. |
| @@ -115,8 +173,11 @@ function _log1pmx_ker(x::Float64) | |||
| end | |||
|
|
|||
|
|
|||
| ## logsumexp | |||
| """ | |||
| logsumexp(x::Real, y::Real) | |||
There was a problem hiding this comment.
we should probably rename this logaddexp?
There was a problem hiding this comment.
Look a bit further down. There is a method for x::AbstractArray{<:Real} for which sum as part of the name is appropriate.
src/misc.jl
Outdated
| @@ -43,7 +72,8 @@ function lstirling_asym(x::Float64) | |||
| end | |||
|
|
|||
| function lstirling_asym(x::Float32) | |||
| t = 1f0/(x*x) | |||
| isinteger(x) && (0 < x ≤ length(lstirlingF32)) && return lstirlingF32[Int(x)] | |||
| t = inv/(abs2(x)) | |||
| @@ -18,3 +18,13 @@ poispdf(λ::Real, x::Real) = exp(poislogpdf(λ, x)) | |||
| poislogpdf(λ::T, x::T) where {T <: Real} = xlogy(x, λ) - λ - lgamma(x + 1) | |||
|
|
|||
| poislogpdf(λ::Number, x::Number) = poislogpdf(promote(float(λ), x)...) | |||
|
|
|||
| #= | |||
| function poislogpdf(λ::Union{Float32,Float64}, x::Union{Float64,Float32,Integer}) | |||
src/misc.jl
Outdated
| function lstirling_asym end | ||
|
|
||
| # For BigInt use the exact calculation of factorial(x-1) | ||
| function lstirling_asym(x::BigInt) |
There was a problem hiding this comment.
I think lgamma is fine everywhere
julia> norm([log(factorial(n)) - lgamma(n + 1) for n in big(1):big(500)])
0.0Maybe better to make change the signature to BigFloat to match the other versions.
There was a problem hiding this comment.
Thanks for catching that.
|
I will add more tests as part of this PR. I just noticed that the author of the |
|
I see now that |
|
regarding the The two most obvious choices seem to be the Lanczos' and Spouge's approximations. They are similar, but the former converges faster (and is used by Boost), whereas the latter can be computed to arbitrary precision (and is used by MPFR). |
|
@simonbyrne My purpose in looking at Do you know of other places where |
|
At the moment, no. But if we can define it for real values we can use it for the gamma and beta densities as well. |
|
I've tweaked the Stirling docs slightly, but lets merge this and we can flesh out other changes later. |
src/misc.jl
Outdated
|
|
||
| lgamma(x) ≈ x*log(x) - x + log(2π/x)/2 = log(x)*(x-0.5) + log2π/2 -x | ||
| ```math | ||
| \log \Gamma(x) \approx x \log(x) - x + log(2π/x)/2 = \log(x)*(x-1/2) + \log(2\pi)/2 - x |
There was a problem hiding this comment.
Missed one log -> \log
My purpose in creating an explicit
BigIntmethod forlstirling_asymwas to get full accuracy on small integer values ofxfrom theFloat64andFloat32methods by evaluatinglog(factorial(x-1))explicitly and storing these in an array. The only downside I can see to this is that the function will be a bit discontinuous at these small integer values (perhaps literally a bit or two of discontinuity).This way the Stirling error function
stirlerr(n)in #33 can be replaced bylstirling_asym(n + 1)