From 520e15411edcd74ae6a9d25148810d0451fdc261 Mon Sep 17 00:00:00 2001 From: Samuel Albert Date: Wed, 21 Apr 2021 14:23:52 -0400 Subject: [PATCH 1/6] added weighted fit for Laplace distribution --- src/univariate/continuous/laplace.jl | 27 +++++++++++++++++++++++---- 1 file changed, 23 insertions(+), 4 deletions(-) diff --git a/src/univariate/continuous/laplace.jl b/src/univariate/continuous/laplace.jl index b47fb0cede..90df543175 100644 --- a/src/univariate/continuous/laplace.jl +++ b/src/univariate/continuous/laplace.jl @@ -114,8 +114,27 @@ rand(rng::AbstractRNG, d::Laplace) = #### Fitting -function fit_mle(::Type{<:Laplace}, x::Array) - xc = copy(x) - a = median!(xc) - Laplace(a, StatsBase.mad!(xc, center=a)) +function fit_mle(::Type{<:Laplace}, x::AbstractArray{T}) where {T <: Real} + μ = median(x) + return Laplace(μ, mean(abs.(x .- μ))) +end + +function fit_mle(::Type{<:Laplace}, x::AbstractArray{T}, w::AbstractArray{T}) where {T <: Real} + sp = sortperm(x) + n = length(x) + sw = sum(w) + highsum = sw + lowsum = zero(T) + idx = 0 + for i = 1:n + lowsum += w[sp[i]] + highsum -= w[sp[i]] + if lowsum >= highsum + idx = sp[i] + break + end + end + μ = x[idx] + θ = sum(w .* abs.(x .- μ)) / sw + return Laplace(μ, θ) end From 6b1da51c86fc0c8b1400207dbd0122faf2c72d00 Mon Sep 17 00:00:00 2001 From: Samuel Albert Date: Wed, 21 Apr 2021 15:34:13 -0400 Subject: [PATCH 2/6] effiency improvements --- src/univariate/continuous/laplace.jl | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/src/univariate/continuous/laplace.jl b/src/univariate/continuous/laplace.jl index 90df543175..bf97f9b2ea 100644 --- a/src/univariate/continuous/laplace.jl +++ b/src/univariate/continuous/laplace.jl @@ -114,8 +114,13 @@ rand(rng::AbstractRNG, d::Laplace) = #### Fitting -function fit_mle(::Type{<:Laplace}, x::AbstractArray{T}) where {T <: Real} +function fit_mle(::Type{<:Laplace}, x::AbstractArray{T}) where {T<:Real} μ = median(x) + θ = zero(T) + for v ∈ x + θ += abs(v - μ) + end + θ /= length(x) return Laplace(μ, mean(abs.(x .- μ))) end @@ -135,6 +140,10 @@ function fit_mle(::Type{<:Laplace}, x::AbstractArray{T}, w::AbstractArray{T}) wh end end μ = x[idx] - θ = sum(w .* abs.(x .- μ)) / sw + θ = zero(T) + for i = 1:length(x) + θ += w[i] * abs(x[i] - μ) + end + θ /= sw return Laplace(μ, θ) end From 059a2377745cf01b9ce605697626b84dcb0f18ae Mon Sep 17 00:00:00 2001 From: Samuel Albert Date: Wed, 21 Apr 2021 15:42:11 -0400 Subject: [PATCH 3/6] bug fix --- src/univariate/continuous/laplace.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/univariate/continuous/laplace.jl b/src/univariate/continuous/laplace.jl index bf97f9b2ea..e853384e0c 100644 --- a/src/univariate/continuous/laplace.jl +++ b/src/univariate/continuous/laplace.jl @@ -121,7 +121,7 @@ function fit_mle(::Type{<:Laplace}, x::AbstractArray{T}) where {T<:Real} θ += abs(v - μ) end θ /= length(x) - return Laplace(μ, mean(abs.(x .- μ))) + return Laplace(μ, θ) end function fit_mle(::Type{<:Laplace}, x::AbstractArray{T}, w::AbstractArray{T}) where {T <: Real} From 857debf0217aae9e38bf0050c7c755cbdac1d970 Mon Sep 17 00:00:00 2001 From: Samuel Albert Date: Wed, 21 Apr 2021 17:07:26 -0400 Subject: [PATCH 4/6] improvements suggested by David Widmann --- src/univariate/continuous/laplace.jl | 38 +++++----------------------- 1 file changed, 6 insertions(+), 32 deletions(-) diff --git a/src/univariate/continuous/laplace.jl b/src/univariate/continuous/laplace.jl index e853384e0c..ebd3945ea8 100644 --- a/src/univariate/continuous/laplace.jl +++ b/src/univariate/continuous/laplace.jl @@ -114,36 +114,10 @@ rand(rng::AbstractRNG, d::Laplace) = #### Fitting -function fit_mle(::Type{<:Laplace}, x::AbstractArray{T}) where {T<:Real} - μ = median(x) - θ = zero(T) - for v ∈ x - θ += abs(v - μ) - end - θ /= length(x) - return Laplace(μ, θ) -end - -function fit_mle(::Type{<:Laplace}, x::AbstractArray{T}, w::AbstractArray{T}) where {T <: Real} - sp = sortperm(x) - n = length(x) - sw = sum(w) - highsum = sw - lowsum = zero(T) - idx = 0 - for i = 1:n - lowsum += w[sp[i]] - highsum -= w[sp[i]] - if lowsum >= highsum - idx = sp[i] - break - end - end - μ = x[idx] - θ = zero(T) - for i = 1:length(x) - θ += w[i] * abs(x[i] - μ) - end - θ /= sw - return Laplace(μ, θ) +function fit_mle(::Type{<:Laplace}, x::AbstractArray{<:Real}) + xc = similar(x) + copyto!(xc, x) + m = median!(xc) + xc .= abs.(x .- m) + return Laplace(m, mean(xc)) end From 0495e6d73d06ec4eb2cffd18c03c2c43434283b2 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Thu, 22 Apr 2021 00:21:53 +0200 Subject: [PATCH 5/6] Lower tolerances in the tests --- test/fit.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/fit.jl b/test/fit.jl index 54be858878..4cf299c6aa 100644 --- a/test/fit.jl +++ b/test/fit.jl @@ -355,10 +355,10 @@ end @testset "Testing fit for Laplace" begin for func in funcs, dist in (Laplace, Laplace{Float64}) - d = fit(dist, func[2](dist(5.0, 3.0), N)) + d = fit(dist, func[2](dist(5.0, 3.0), N + 1)) @test isa(d, dist) - @test isapprox(location(d), 5.0, atol=0.1) - @test isapprox(scale(d) , 3.0, atol=0.2) + @test isapprox(location(d), 5.0, atol=0.02) + @test isapprox(scale(d) , 3.0, atol=0.02) end end From e3c70c001800f37b6b47439d603eecc041b4808a Mon Sep 17 00:00:00 2001 From: David Widmann Date: Thu, 22 Apr 2021 00:22:11 +0200 Subject: [PATCH 6/6] Bump version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 8f3b50ac91..32ec76710a 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Distributions" uuid = "31c24e10-a181-5473-b8eb-7969acd0382f" authors = ["JuliaStats"] -version = "0.24.17" +version = "0.24.18" [deps] FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"