From 09e995509ef60a62c87c99423da6db79f17de42f Mon Sep 17 00:00:00 2001 From: pagnani Date: Wed, 6 Feb 2019 11:40:14 -0500 Subject: [PATCH] added dependency from Distributions. added a max to set to zero negative variances: issue #827 in Distributions.jl --- src/priors.jl | 33 ++++++--------------------------- 1 file changed, 6 insertions(+), 27 deletions(-) diff --git a/src/priors.jl b/src/priors.jl index f288abe..e0d4c3a 100644 --- a/src/priors.jl +++ b/src/priors.jl @@ -1,4 +1,5 @@ using FastGaussQuadrature, ForwardDiff +import Distributions: TruncatedNormal, mean, var Φ(x) = 0.5*(1+erf(x/sqrt(2.0))) ϕ(x) = exp(-x.^2/2)/sqrt(2π) @@ -45,35 +46,13 @@ struct IntervalPrior{T<:Real} <: Prior u::T end -function moments(p0::IntervalPrior,μ,σ) - xl = (p0.l - μ)/σ - xu = (p0.u - μ)/σ - minval = min(abs(xl), abs(xu)) - - if xu - xl < 1e-10 - return 0.5 * (xu + xl), -1 - end - - if minval <= 6.0 || xl * xu <= 0 - ϕu, Φu, ϕl, Φl = ϕ(xu), Φ(xu), ϕ(xl), Φ(xl) - av = (ϕl - ϕu) / (Φu - Φl) - mom2 = (xl * ϕl - xu * ϕu) / (Φu - Φl) - else - Δ = (xu^2 - xl^2) * 0.5 - if Δ > 40.0 - av = xl^5 / (3.0 - xl^2 + xl^4) - mom2 = xl^6 / (3.0 - xl^2 + xl^4) - else - eΔ = exp(Δ) - av = (xl * xu)^5 * (1. - eΔ) / (-eΔ * (3.0 - xl^2 + xl^4) * xu^5 + xl^5 * (3.0 - xu^2 + xu^4)) - mom2 = (xl * xu)^5 * (xu - xl * eΔ) / (-eΔ * (3.0 - xl^2 + xl^4) * xu^5 + xl^5 * (3.0 - xu^2 + xu^4)) - end - end - va = mom2 - av^2 - return μ + av * σ, σ^2 * (1 + va) +function moments(p0::IntervalPrior{T},μ,σ) where T<:Real + pr = TruncatedNormal(μ,σ,p0.l,p0.u) + return mean(pr),max(zero(T),var(pr)) + # the max fix an error in Distributions + # issue #827 in Distributions.jl end - """ Spike-and-slab prior