Repeated draws from Multinomial causes: ERROR: prob must be in [0, 1]

[lejon]

Hi,

I think that I might have found a bug in the multinomial implementation. It seems that sometimes the binomial is given a p > 1 which I believe is caused by a rounding error. The problem can be reproduced with the below code:

$cat distr_bug.jl

using Distributions
alphas = fill(0.1, 10)
dir = Dirichlet(alphas)

cnt = 0
while(true) 
    dirdraw = rand(dir)
    if( mod(cnt, 1000) == 0)
        println("$cnt DIR: $dirdraw")
    end
    mul = Multinomial(1,dirdraw)
    muldraw = rand(mul)
    cnt += 1
end

$ julia distr_bug.jl

ERROR: prob must be in [0, 1]
 in Binomial at /Users/lejon/.julia/Distributions/src/univariate/binomial.jl:7
 in multinom_rand! at /Users/lejon/.julia/Distributions/src/multivariate/multinomial.jl:135
 in anonymous at no file:12
 in include at boot.jl:238
 in include_from_node1 at loading.jl:119
 in process_options at client.jl:307
 in _start at client.jl:387
at /Users/lejon/Data/coding/julia/distr_bug.jl:14

-Leif

Edit (kms): added formatting with triple ticks, for code readability

[kmsquire]

I tried to fix this the hard way, borrowing KBN summation from cumsum_kbn, but that didn't seem to fix it:

function multinom_rand!{T<:Real}(n::Int, p::Vector{Float64}, x::AbstractVector{T})
    ...
    while i < km1 && n > 0
        i += 1
        @inbounds pi = p[i]
        xi = rand(Binomial(n, pi / (rp+c)))
        @inbounds x[i] = xi
        n -= xi
        rp1 = rp - pi
        if abs(rp) >= abs(pi)
            c += ((rp-rp1) - pi)
        else
            c += (rp - (rp1+pi))
        end
        rp = rp1
    end
    ...
end

julia> cnt = 0; while(true) 
          dirdraw = rand(dir)
          if( mod(cnt, 1000) == 0)
             println("$cnt DIR: $dirdraw")
          end
          mul = Multinomial(1,dirdraw)
          muldraw = rand(mul)
          cnt += 1
       end
0 DIR: 8656205104069296e-27
.0018630726231783187
.0034847262841801764
.07958155215732116
.0048209869575798914
.492338714296625
.29905579505252033
11247520370455625e-23
.0031099475237326697
.1157450926210026

ERROR: prob must be in [0, 1] (got p=1.0000000000000002)
 in error at error.jl:21
 in Binomial at /home/kmsquire/.julia/v0.3/Distributions/src/univariate/binomial.jl:6
 in multinom_rand! at /home/kmsquire/.julia/v0.3/Distributions/src/multivariate/multinomial.jl:136
 in anonymous at no file:7

So, it looks like something more direct is needed (e.g., Binomial(n, min(pi / rp, 1.0))))

This still might cause problems if rp is ever less than zero (is that possible?), or alternatively, if rp > sum(p) which always causes probability mass to collect in the last position (which is extremely unlikely to actually cause any problems, but is theoretically possible.)

[johnmyleswhite]

Let's use the Binomial(n, min(pi / rp, 1.0)) fix for now. I don't have enough feel for floating point arithmetic to make this more stable otherwise. Maybe Dahua will have some insight?

[johnmyleswhite]

64625eccf1186ddda1615bfeeb749c56bd641dbf offers a quick fix. It would be good to see if we can come up with something cleaner, but this change seems almost uniformly better than just failing.

[johnmyleswhite]

Anyone have more thoughts on this? It worries me that we're getting floating point errors that break things, but I don't see an easy solution to this one. Maybe @StefanKarpinski will have some ideas.

[lindahua]

I am looking into this.

[lindahua]

I found a case that can reliably reproduce this issue:

p = [
4182606124264053e-42
.015712938543046305
3600131090855299e-20
36742804055163726e-21
.3052582444314661
8005659074715409e-22
.06627521733229433
.025244201528410688
.5874358534839115
11244228473164858e-52]

Here is how the function proceeds:

i = 1, pi = 4.182606124264053e-27, rp = 1.0
i = 2, pi = 0.015712938543046305, rp = 1.0
i = 3, pi = 3.600131090855299e-5, rp = 0.9842870614569537
i = 4, pi = 3.6742804055163726e-5, rp = 0.9842510601460451
i = 5, pi = 0.3052582444314661, rp = 0.9842143173419899
i = 6, pi = 8.005659074715409e-7, rp = 0.6789560729105238
i = 7, pi = 0.06627521733229433, rp = 0.6789552723446164
i = 8, pi = 0.025244201528410688, rp = 0.612680055012322
i = 9, pi = 0.5874358534839115, rp = 0.5874358534839114

That's how pi ends up greater than rp (by a extremely small margin).

I am working on a better implementation of the sampling procedure.

[johnmyleswhite]

I'm hesitant to do this, but I wondered yesterday while making our quick fix if perhaps Binomial should tolerate values of p that are within a small epsilon on 0 and 1, transforming them into 0 and 1 respectively.

[lindahua]

This issue can be circumvented without changing the behavior of Binomial.

[lindahua]

A revised implementation in bd9cb6c.

Please take a look at it, and see if we need further efforts. We may close the issue if this version is satisfactory.

[lindahua]

I made some cosmetic modification. I think the current implementation is good enough for most purposes. I am closing this issue. Please feel free to reopen it if any problems come up later.

[StefanKarpinski]

Yeah, I don't know. Maxing out at 1.0 is what I've done when I've encountered this in the past. Not sure that it's the best solution possible, but it's certainly better than failing.