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

Distributions.jl/src/samplers/vonmisesfisher.jl

Line 103 in 8887d7a

v[1] /= s

: 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

Distributions.jl/src/samplers/vonmisesfisher.jl

Line 103 in 8887d7a

v[1] /= s

: 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

Distributions.jl/src/samplers/vonmisesfisher.jl

Lines 87 to 102 in 13029c0

	function _vmf_householder_vec(μ::Vector{Float64})
	# assuming μ is a unit-vector (which it should be)
	# can compute v in a single pass over μ
	p = length(μ)
	v = similar(μ)
	v[1] = μ[1] - 1.0
	s = sqrt(-2*v[1])
	v[1] /= s
	@inbounds for i in 2:p
	v[i] = μ[i] / s
	end
	return v
	end

for the case that s is numerically close to zero. That should give the correct (identity) rotation in

Distributions.jl/src/samplers/vonmisesfisher.jl

Line 46 in 13029c0

return _vmf_rot!(spl.v, x)

Maybe @emerali has a suggestion?