Add Landau distribution #1245
Add Landau distribution #1245
Conversation
Codecov Report
@@ Coverage Diff @@
## master #1245 +/- ##
==========================================
- Coverage 81.83% 80.66% -1.17%
==========================================
Files 117 118 +1
Lines 6589 6689 +100
==========================================
+ Hits 5392 5396 +4
- Misses 1197 1293 +96
Continue to review full report at Codecov.
|
| exp(im*t*d.μ - 2im*d.θ*t/π*log(abs(t)) - d.θ*abs(t)) | ||
| end | ||
|
|
||
| # https://root.cern.ch/doc/v622/PdfFuncMathCore_8cxx_source.html |
There was a problem hiding this comment.
I'm not an expert but I think this is problematic since root uses the LGPL license: https://root.cern/about/license/
There was a problem hiding this comment.
yes, someone pointed out on Slack as well, it is true we cannot put derived code directly. so this may call for a new pkg
| Landau(μ::T) where {T <: Real} = Landau(μ, one(T)) | ||
| Landau() = Landau(0.0, 1.0, check_args=false) |
There was a problem hiding this comment.
| Landau(μ::T) where {T <: Real} = Landau(μ, one(T)) | |
| Landau() = Landau(0.0, 1.0, check_args=false) | |
| Landau(μ::Real = 0) = Landau(μ, one(μ); check_args=false) |
| @@ -0,0 +1,191 @@ | |||
| using StaticArrays | |||
There was a problem hiding this comment.
I assume dependencies should be loaded in src/Distributions.jl. I'm a bit surprised though that Distributions depends on StaticArrays but StaticArrays doesn't seem to be imported.
There was a problem hiding this comment.
I think it only uses them in test, so technically we should change that in Project.toml to reflect this
| function cf(d::Landau, t::Real) | ||
| exp(im*t*d.μ - 2im*d.θ*t/π*log(abs(t)) - d.θ*abs(t)) | ||
| end |
There was a problem hiding this comment.
The following definition uses cis(x) = exp(im*x) and avoids NaN values for t = 0:
| function cf(d::Landau, t::Real) | |
| exp(im*t*d.μ - 2im*d.θ*t/π*log(abs(t)) - d.θ*abs(t)) | |
| end | |
| function cf(d::Landau, t::Real) | |
| z = t * (d.μ + d.θ * (sign(t) * im - twoinvπ * log(abs(t)))) | |
| return iszero(t) ? cis(zero(z)) : cis(z) | |
| end |
| # https://root.cern.ch/doc/v622/PdfFuncMathCore_8cxx_source.html | ||
| function _landau_pdf(x, xi=1, x0=0) | ||
| # landau pdf : algorithm from cernlib g110 denlan | ||
| p1 = @SVector [0.4259894875,-0.1249762550, 0.03984243700, -0.006298287635, 0.001511162253] |
There was a problem hiding this comment.
In general, I think it might be easier to just use evalpoly or @evalpoly.
|
|
||
| if (v < -5.5) | ||
| u = exp(v+1) | ||
| lan = 0.3989422803*exp(-1. /u)*sqrt(u)*(1+(a1[2]+(a1[3]+a1[4]*u)*u)*u) |
There was a problem hiding this comment.
Same here, I guess using evalpoly or @evalpoly would be easier to read.
Co-authored-by: David Widmann <devmotion@users.noreply.github.com>
Co-authored-by: David Widmann <devmotion@users.noreply.github.com>
Co-authored-by: David Widmann <devmotion@users.noreply.github.com>
3 participants
I don't know what's the standard for numerical methods, I'd happy to split into a separate package if this is not suitable for Distributions.jl