mode(Binomial) (docs?) wrong when (n+1)*p is integer

[andreasnoack]

When (n+1)*p is integer, then there are two modes, see https://onlinelibrary.wiley.com/doi/10.1111/j.1467-9574.1980.tb00681.x. The documentation states

  mode(d::UnivariateDistribution)

  Returns the first mode.

but the second is returned, e.g.

julia> [(i, (pdf(Binomial(5, 4/6), i))) for i in 0:5]
6-element Vector{Tuple{Int64, Float64}}:
 (0, 0.004115226337448564)
 (1, 0.04115226337448562)
 (2, 0.16460905349794255)
 (3, 0.329218106995885)
 (4, 0.3292181069958847)
 (5, 0.13168724279835395)

julia> mode(Binomial(5, 4/6))
4

[marcusps]

Also, modes(d) returns [4] in this case, when it should be [3, 4]