Prevent rejection sampling if truncated normal variance is 0

[quildtide]

Fix #1264

Upon initially investigating this, I discovered that this was not a crash, but instead a process that either does not resolve or a process that takes an extremely long time to resolve. My initial suspicion was that this was an infinite loop. Furthermore, it can be triggered with rand(truncate(Normal(0,0), x, y)) where x and y are any numbers that are 0 or of different sign, along with some other cases that I won't immediately delve into. I suspect that the mean of the Normal distribution is also irrelevant, but I did not bother testing this.

We have a method for rand(d::Truncated{Normal{T}, Continus}) defined at line 140 of truncated/normal.jl (

Distributions.jl/src/truncated/normal.jl

Line 140 in 863844c

function rand(rng::AbstractRNG, d::Truncated{Normal{T},Continuous}) where T <: Real

):

function rand(rng::AbstractRNG, d::Truncated{Normal{T},Continuous}) where T <: Real
    d0 = d.untruncated
    μ = mean(d0)
    σ = std(d0)
    if isfinite(μ)
        a = (d.lower - μ) / σ
        b = (d.upper - μ) / σ
        z = randnt(rng, a, b, d.tp)
        return μ + σ * z
    else
        return clamp(μ, d.lower, d.upper)
    end
end

isfinite(mean(truncated(Normal(0,0), 0, 1))) correctly evaluates to true, bringing us to the randnt function (although it's worth noting that we're passing exclusively NaN, Inf, and -Inf to the randnt function. This is because we are dividing by the standard deviation, which is 0, in order to standardize the bounds.

Within the 2nd part of randnt, the first two cases, when triggered, yield acceptable results even for Normal(0,0), e.g. they always yield 0.0 within a reasonable amount of time.

When neither of these cases are triggered, a third case is triggered which, under a wide variety of x,y where our underlying distribution is Normal(0,0) and our bounds are x,y, never exits. Under these circumstances, rho is always NaN and we can only exit the loop when r < rho (where r is randomly generated). r < NaN will never evaluate to true, so we never exit the loop.

---

However, there's no need to actually delve deep into randnt to fix this. We should probably just avoid trying to feed randnt the "standardized" values from a 0-variance distribution in the first place.

I propose fixing this at the level of rand(rng::AbstractRNG, d::Truncated{Normal{T},Continuous}) where T <: Real:

We return the mean when a user attempts to sample from a 0-variance truncated normal distribution.

[quildtide]

On the other hand, an issue with my "fix":

julia> rand(Distributions.truncated(Distributions.Normal(0, 0), 0, Inf))
0.0

julia> rand(Distributions.truncated(Distributions.Normal(-1, 0), 0, Inf))
-1.0

It's probably also necessary to prevent (via error) the creation of something like Distributions.truncated(Distributions.Normal(-1, 0), 0, Inf) in the first place, since the distribution has no support at all. At the least, there should be a warning.

On the other hand, this behavior is consistent with other behavior currently observed with truncated distributions with no remaining support, such as:

julia> rand(Distributions.truncated(Distributions.Beta(1, 1), 2, 3))
1.0

julia> rand(Distributions.truncated(Distributions.Beta(1, 1), -5, -1))
0.0

julia> rand(Distributions.truncated(Distributions.Binomial(1, 1), -5, -1))
0

julia> rand(Distributions.truncated(Distributions.Binomial(1, 1), 2, 3))
1

This has been mentioned before back in 2019 in #843, but it does not seem like anything has been done since then.