Commit 5ae0432c authored by Dahua Lin's avatar Dahua Lin

new MLE estimation for von-Mises Fisher distributions

parent 464e8701
......@@ -65,3 +65,47 @@ _rand!(d::VonMisesFisher, x::DenseVector) = _rand!(sampler(d), x)
_rand!(d::VonMisesFisher, x::DenseMatrix) = _rand!(sampler(d), x)
### Estimation
function fit_mle(::Type{VonMisesFisher}, X::Matrix{Float64})
r = vec(sum(X, 2))
n = size(X, 2)
r_nrm = vecnorm(r)
μ = scale!(r, 1.0 / r_nrm)
ρ = r_nrm / n
κ = _vmf_estkappa(length(μ), ρ)
VonMisesFisher(μ, κ)
end
fit_mle{T<:Real}(::Type{VonMisesFisher}, X::Matrix{T}) = fit_mle(VonMisesFisher, float64(X))
function _vmf_estkappa(p::Int, ρ::Float64)
# Using the fixed-point iteration algorithm in the following paper:
#
# Akihiro Tanabe, Kenji Fukumizu, and Shigeyuki Oba, Takashi Takenouchi, and Shin Ishii
# Parameter estimation for von Mises-Fisher distributions.
# Computational Statistics, 2007, Vol. 22:145-157.
#
const maxiter = 200
half_p = 0.5 * p
ρ2 = abs2(ρ)
κ = ρ * (p - ρ2) / (1 - ρ2)
i = 0
while i < maxiter
i += 1
κ_prev = κ
a = (ρ / _vmfA(half_p, κ))
# println("i = $i, a = $a, abs(a - 1) = $(abs(a - 1))")
κ *= a
if abs(a - 1.0) < 1.0e-12
break
end
end
return κ
end
_vmfA(half_p::Float64, κ::Float64) = besseli(half_p, κ) / besseli(half_p - 1.0, κ)
......@@ -15,7 +15,7 @@ function gen_vmf_tdata(n::Int, p::Int)
return X
end
function test_vonmisesfisher(p::Int, κ::Float64, n::Int)
function test_vonmisesfisher(p::Int, κ::Float64, n::Int, ns::Int)
μ = randn(p)
μ = μ ./ vecnorm(μ)
......@@ -50,19 +50,26 @@ function test_vonmisesfisher(p::Int, κ::Float64, n::Int)
for i = 1:n
@test_approx_eq vecnorm(X[:,i]) 1.0
end
# MLE
X = rand(d, ns)
d_est = fit_mle(VonMisesFisher, X)
@test isa(d_est, VonMisesFisher)
@test_approx_eq_eps d_est.μ μ 0.01
@test_approx_eq_eps d_est.κ κ κ * 0.01
end
## General testing
n = 1000
ns = 10^6
for (p, κ) in [(2, 1.0),
(2, 5.0),
(3, 1.0),
(3, 5.0),
(5, 2.0)]
test_vonmisesfisher(p, κ, n)
test_vonmisesfisher(p, κ, n, ns)
end
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment