Merge pull request #223 from jdrugo/addvonmises

Add support for univariate von Mises distribution

doc/source/univariate.rst

@@ -219,6 +219,7 @@ List of Distributions
     - :ref:`rayleigh`
     - :ref:`tdist`
     - :ref:`uniform`
+    - :ref:`vonmises`
     - :ref:`weibull`
 
 
@@ -745,6 +746,24 @@ The probability density function of a `Continuous Uniform distribution <http://e
     Uniform()        # Uniform distribution over [0, 1]
     Uniform(a, b)    # Uniform distribution over [a, b]
 
+.. _vonmises:
+
+Von Mises Distribution
+~~~~~~~~~~~~~~~~~~~~~~~
+
+The probability density function of a `von Mises distribution <http://en.wikipedia.org/wiki/Von_Mises_distribution>`_ with mean μ and concentration κ is
+
+.. math::
+
+    f(x; \mu, \kappa) = \frac{1}{2 \pi I_0(\kappa)} \exp \left( \kappa \cos (x - \mu) \right)
+
+.. code-block:: julia
+
+    VonMises()       # von Mises distribution with zero mean and unit concentration
+    VonMises(κ)      # von Mises distribution with zero mean and concentration κ
+    VonMises(μ, κ)   # von Mises distribution with mean μ and concentration κ
+
+
 .. _weibull:
 
 Weibull Distribution

src/Distributions.jl

@@ -105,6 +105,7 @@ export
     SymTriangularDist,
     Truncated,
     Uniform,
+    VonMises,
     VonMisesFisher,
     Weibull,
     Wishart,
@@ -117,6 +118,10 @@ export
     cdf,                # cumulative distribution function
     cf,                 # characteristic function
     cgf,                # cumulant generating function
+    circmean,           # mean of circular distribution
+    circmedian,         # median of circular distribution
+    circmode,           # mode of circular distribution
+    circvar,            # variance of circular distribution
     cquantile,          # complementary quantile (i.e. using prob in right hand tail)
     cumulant,           # cumulants of distribution
     complete,           # turn an incomplete formulation into a complete distribution

src/samplers.jl

@@ -5,7 +5,8 @@ for fname in ["categorical.jl",
               "poisson.jl", 
               "exponential.jl",
               "gamma.jl", 
-              "multinomial.jl"]
+              "multinomial.jl",
+              "vonmises.jl"]
 
     include(joinpath("samplers", fname))
 end

src/samplers/vonmises.jl

@@ -0,0 +1,50 @@
+
+immutable VonMisesSampler <: Sampleable{Univariate,Continuous}
+    μ::Float64
+    κ::Float64
+    r::Float64
+end
+
+# algorithm from
+#     DJ Best & NI Fisher (1979). Efficient Simulation of the von Mises
+#     Distribution. Journal of the Royal Statistical Society. Series C
+#     (Applied Statistics), 28(2), 152-157.
+function rand(s::VonMisesSampler)
+    f = 0.0
+    while true
+        t, u = 0.0, 0.0
+        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)
+                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
+    end
+    mod(s.μ + (rand() > 0.5 ? acos(f) : -acos(f)), twoπ)
+end
+
+immutable VonMisesNormalApproxSampler <: Sampleable{Univariate,Continuous}
+    μ::Float64
+    σ::Float64
+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
+

src/univariate/vonmises.jl

@@ -0,0 +1,62 @@
+# von Mises distribution
+
+immutable VonMises <: ContinuousUnivariateDistribution
+    μ::Float64  # mean
+    κ::Float64  # concentration
+
+    function VonMises(μ::Real, κ::Real)
+        κ > zero(κ) || error("kappa must be positive")
+        new(float64(μ), float64(κ))
+    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.κ)
+
+function entropy(d::VonMises)
+	I0κ = besselix(0.0, d.κ)
+	log(twoπ * I0κ) - d.κ * (besselix(1, d.κ) / I0κ - 1.0)
+end
+
+## 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, κ))
+end
+
+## Sampling
+sampler(d::VonMises) = VonMisesSampler(d.μ, d.κ)

src/univariates.jl

@@ -153,6 +153,7 @@ for finame in ["arcsine.jl",
                "tdist.jl",
                "symtriangular.jl",
                "uniform.jl",
+               "vonmises.jl",
                "weibull.jl"
                ]
     include(joinpath("univariate", finame))