Ce serveur Gitlab sera éteint le 30 juin 2020, pensez à migrer vos projets vers les serveurs gitlab-research.centralesupelec.fr et gitlab-student.centralesupelec.fr !

Commit b830ba9d authored by Dahua Lin's avatar Dahua Lin

Merge pull request #300 from JuliaStats/dh/vonmise2

Reimplement VonMise distribution. Fix #273
parents 8dc6124c d2bfb7da
......@@ -3,6 +3,12 @@ immutable VonMisesSampler <: Sampleable{Univariate,Continuous}
μ::Float64
κ::Float64
r::Float64
function VonMisesSampler(μ::Float64, κ::Float64)
τ = 1.0 + sqrt(1.0 + 4 * abs2(κ))
ρ = (τ - sqrt(2.0 * τ)) / (2.0 * κ)
new(μ, κ, (1.0 + abs2(ρ)) / (2.0 * ρ))
end
end
# algorithm from
......@@ -11,40 +17,30 @@ end
# (Applied Statistics), 28(2), 152-157.
function rand(s::VonMisesSampler)
f = 0.0
while true
t, u = 0.0, 0.0
x::Float64
if s.κ > 700.0
x = s.μ + randn() / sqrt(s.κ)
else
while true
const v, w = rand() - 0.5, rand() - 0.5
const d, e = v ^ 2, w ^ 2
if d + e <= 0.25
t = d / e
u = 4 * (d + e)
t, u = 0.0, 0.0
while true
d = abs2(rand() - 0.5)
e = abs2(rand() - 0.5)
if d + e <= 0.25
t = d / e
u = 4 * (d + e)
break
end
end
z = (1.0 - t) / (1.0 + t)
f = (1.0 + s.r * z) / (s.r + z)
c = s.κ * (s.r - f)
if c * (2.0 - c) > u || log(c / u) + 1 >= c
break
end
end
const z = (1.0 - t) / (1.0 + t)
f = (1.0 + s.r * z) / (s.r + z)
const c = s.κ * (s.r - f)
if c * (2.0 - c) > u || log(c / u) + 1 >= c
break
end
acf = acos(f)
x = s.μ + (randbool() ? acf : -acf)
end
mod(s.μ + (rand() > 0.5 ? acos(f) : -acos(f)), twoπ)
end
immutable VonMisesNormalApproxSampler <: Sampleable{Univariate,Continuous}
μ::Float64
σ::Float64
return x
end
rand(s::VonMisesNormalApproxSampler) = mod(s.μ + s.σ * randn(), twoπ)
# normal approximation for large concentrations
VonMisesSampler(μ::Float64, κ::Float64) =
κ > 700.0 ? VonMisesNormalApproxSampler(μ, sqrt(1.0 / κ)) :
begin
τ = 1.0 + sqrt(1.0 + 4 * κ ^ 2)
ρ = (τ - sqrt(2.0 * τ)) / (2.0 * κ)
VonMisesSampler(μ, κ, (1.0 + ρ ^ 2) / (2.0 * ρ))
end
......@@ -155,8 +155,15 @@ function test_samples(s::Sampleable{Univariate, Continuous}, # the sampleable
vmin = minimum(distr)
vmax = maximum(distr)
rmin = quantile(distr, 0.01)
rmax = quantile(distr, 0.99)
rmin::Float64
rmax::Float64
if applicable(quantile, distr, 0.5)
rmin = quantile(distr, 0.01)
rmax = quantile(distr, 0.99)
elseif isfinite(vmin) && isfinite(vmax)
rmin = vmin
rmax = vmax
end
edges = linspace(rmin, rmax, nbins+1)
# determine confidence intervals for counts:
......
# von Mises distribution
#
# Implemented based on Wikipedia
#
immutable VonMises <: ContinuousUnivariateDistribution
μ::Float64 # mean
κ::Float64 # concentration
μ::Float64 # mean
κ::Float64 # concentration
I0κ::Float64 # I0(κ), where I0 is the modified Bessel function of order 0
function VonMises(μ::Real, κ::Real)
κ > zero(κ) || error("kappa must be positive")
new(float64(μ), float64(κ))
new(float64(μ), float64(κ), besseli(0, κ))
end
VonMises(κ::Real) = VonMises(0.0, float64(κ))
VonMises() = VonMises(0.0, 1.0)
end
## Properties
circmean(d::VonMises) = d.μ
circmedian(d::VonMises) = d.μ
circmode(d::VonMises) = d.μ
circvar(d::VonMises) = 1.0 - besselix(1, d.κ) / besselix(0, d.κ)
show(io::IO, d::VonMises) = show(io, d, (:μ, :κ))
function entropy(d::VonMises)
I0κ = besselix(0.0, d.κ)
log(twoπ * I0κ) - d.κ * (besselix(1, d.κ) / I0κ - 1.0)
end
### Properties
minimum(d::VonMises) = d.μ - π
maximum(d::VonMises) = d.μ + π
islowerbounded(d::VonMises) = true
isupperbounded(d::VonMises) = true
insupport(d::VonMises, x::Real) = -π <= (x - d.μ) <= π
mean(d::VonMises) = d.μ
median(d::VonMises) = d.μ
mode(d::VonMises) = d.μ
circvar(d::VonMises) = 1.0 - besseli(1, d.κ) / I0κ
entropy(d::VonMises) = log(twoπ * d.I0κ) - d.κ * (besseli(1, d.κ) / d.I0κ)
cf(d::VonMises, t::Real) = (besseli(abs(t), d.k) / I0κ) * exp(im * t * d.μ)
## Functions
pdf(d::VonMises, x::Real) = exp(d.κ * (cos(x - d.μ) - 1.0)) / (twoπ * besselix(0, d.κ))
logpdf(d::VonMises, x::Real) = d.κ * (cos(x - d.μ) - 1.0) - log2π - log(besselix(0, d.κ))
cf(d::VonMises, t::Real) = besselix(abs(t), d.k) / besselix(0.0, d.κ) * exp(im * t * d.μ)
cdf(d::VonMises, x::Real) = cdf(d, x, d.μ - π)
function cdf(d::VonMises, x::Real, from::Real)
const tol = 1.0e-20
x = mod(x - from, twoπ)
μ = mod(d.μ - from, twoπ)
if μ == 0.0
return vmcdfseries(d.κ, x, tol)
elseif x <= μ
upper = mod(x - μ, twoπ)
if upper == 0.0
upper = twoπ
end
return vmcdfseries(d.κ, upper, tol) - vmcdfseries(d.κ, mod(-μ, twoπ), tol)
else
return vmcdfseries(d.κ, x - μ, tol) - vmcdfseries(d.κ, μ, tol)
end
end
function vmcdfseries(κ::Real, x::Real, tol::Real)
j, s = 1, 0.0
while true
sj = besselix(j, κ) * sin(j * x) / j
s += sj
j += 1
abs(sj) >= tol || break
end
x / twoπ + s / (π * besselix(0, κ))
### Functions
pdf(d::VonMises, x::Real) = exp(d.κ * cos(float64(x) - d.μ)) / (twoπ * d.I0κ)
logpdf(d::VonMises, x::Real) = d.κ * cos(float64(x) - d.μ) - log(d.I0κ) - log2π
cdf(d::VonMises, x::Real) = _vmcdf(d.κ, d.I0κ, float64(x) - d.μ, 1.0e-15)
function _vmcdf(κ::Float64, I0κ::Float64, x::Float64, tol::Float64)
j = 1
cj = besseli(j, κ) / j
s = cj * sin(j * x)
while abs(cj) > tol
j += 1
cj = besseli(j, κ) / j
s += cj * sin(j * x)
end
return (x + 2.0 * s / I0κ) / twoπ + 0.5
end
## Sampling
......
......@@ -108,5 +108,15 @@ for distr in [
println(" testing $(distr)")
test_distr(distr, n_tsamples)
end
end
for distr in [
VonMises(0.0, 1.0),
VonMises(0.5, 1.0),
VonMises(0.5, 2.0) ]
println(" testing $(distr)")
test_samples(distr, n_tsamples)
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