Sampling from von Mises Fisher yields NaN

[lanceXwq]

This happens when the mean direction unit vector μ is along the x-axis. For instance:

julia> vMF = VonMisesFisher([1, 0], 1.0)
VonMisesFisher{Float64}(μ=[1.0, 0.0], κ=1.0)
julia> rand(vMF)
2-element Vector{Float64}:
 NaN
 NaN

This problem seems to persist in higher dimensions.

[devmotion]

I'm not familiar with the implementation of the VonMisesFisher sampler but I suspect that the NaNs arise in https://github.com/JuliaStats/Distributions.jl/blob/8887d7a16a2801dc869ab42ba03cad03efc5ed9d/src/samplers/vonmisesfisher.jl#L103: for μ=[1.0, 0.0] we have v[1] = s = 0 in this line, and hence v[1] /= s should result in v[1] = NaN.

[trahflow]

I just stumbled over this as well.

I'm not familiar with the implementation of the VonMisesFisher sampler but I suspect that the NaNs arise in

https://github.com/JuliaStats/Distributions.jl/blob/8887d7a16a2801dc869ab42ba03cad03efc5ed9d/src/samplers/vonmisesfisher.jl#L103 : for μ=[1.0, 0.0] we have v[1] = s = 0 in this line, and hence v[1] /= s should result in v[1] = NaN.

Yes, this is what happens. It seems this was introduced in #1162. At least a quick check shows the previous version which constructs an explicit rotation matrix does not suffer from this.

A very quick, albeit a bit dirty fix would be to return a zero vector in

https://github.com/JuliaStats/Distributions.jl/blob/13029c03fa885a2b4512b954e61f9d5a7dfa0612/src/samplers/vonmisesfisher.jl#L87-L102 for the case that s is numerically close to zero. That should give the correct (identity) rotation in https://github.com/JuliaStats/Distributions.jl/blob/13029c03fa885a2b4512b954e61f9d5a7dfa0612/src/samplers/vonmisesfisher.jl#L46 Maybe @emerali has a suggestion?