    {-# LANGUAGE BangPatterns, FlexibleContexts, UnboxedTuples #-}
    -- |
    -- Module    : Statistics.Sample.KernelDensity
    -- Copyright : (c) 2011 Bryan O'Sullivan
    -- License   : BSD3
    --
    -- Maintainer  : bos@serpentine.com
    -- Stability   : experimental
    -- Portability : portable
    --
    -- Kernel density estimation.  This module provides a fast, robust,
    -- non-parametric way to estimate the probability density function of
    -- a sample.
    --
    -- This estimator does not use the commonly employed \"Gaussian rule
    -- of thumb\".  As a result, it outperforms many plug-in methods on
    -- multimodal samples with widely separated modes.

    module Statistics.Sample.KernelDensity
        (
        -- * Estimation functions
          kde
        , kde_
        -- * References
        -- $references
        ) where

    import Data.Default.Class
    import Numeric.MathFunctions.Constants (m_sqrt_2_pi)
    import Numeric.RootFinding             (fromRoot, ridders, RiddersParam(..), Tolerance(..))
    import Prelude hiding (const, min, max, sum)
    import Statistics.Function (minMax, nextHighestPowerOfTwo)
    import Statistics.Sample.Histogram (histogram_)
    import Statistics.Sample.Internal (sum)
    import Statistics.Transform (CD, dct, idct)
    import qualified Data.Vector.Generic  as G
    import qualified Data.Vector.Unboxed  as U
    import qualified Data.Vector          as V

    -- | Gaussian kernel density estimator for one-dimensional data, using
    -- the method of Botev et al.
    --
    -- The result is a pair of vectors, containing:
    --
    -- * The coordinates of each mesh point.  The mesh interval is chosen
    --   to be 20% larger than the range of the sample.  (To specify the
    --   mesh interval, use 'kde_'.)
    --
    -- * Density estimates at each mesh point.
    kde :: (G.Vector v CD, G.Vector v Double, G.Vector v Int)
        => Int
        -- ^ The number of mesh points to use in the uniform discretization
        -- of the interval @(min,max)@.  If this value is not a power of
        -- two, then it is rounded up to the next power of two.
        -> v Double -> (v Double, v Double)
    kde :: forall (v :: * -> *).
    (Vector v CD, Vector v Double, Vector v Int) =>
    Int -> v Double -> (v Double, v Double)
    kde Int
    n0 v Double
    xs = Int -> Double -> Double -> v Double -> (v Double, v Double)
    forall (v :: * -> *).
    (Vector v CD, Vector v Double, Vector v Int) =>
    Int -> Double -> Double -> v Double -> (v Double, v Double)
    kde_ Int
    n0 (Double
    lo Double -> Double -> Double
    forall a. Num a => a -> a -> a
    - Double
    range Double -> Double -> Double
    forall a. Fractional a => a -> a -> a
    / Double
    10) (Double
    hi Double -> Double -> Double
    forall a. Num a => a -> a -> a
    + Double
    range Double -> Double -> Double
    forall a. Fractional a => a -> a -> a
    / Double
    10) v Double
    xs
      where
        (Double
    lo,Double
    hi) = v Double -> (Double, Double)
    forall (v :: * -> *).
    Vector v Double =>
    v Double -> (Double, Double)
    minMax v Double
    xs
        range :: Double
    range   | v Double -> Int
    forall (v :: * -> *) a. Vector v a => v a -> Int
    G.length v Double
    xs Int -> Int -> Bool
    forall a. Ord a => a -> a -> Bool
    <= Int
    1 = Double
    1       -- Unreasonable guess
                | Double
    lo Double -> Double -> Bool
    forall a. Eq a => a -> a -> Bool
    == Double
    hi         = Double
    1       -- All elements are equal
                | Bool
    otherwise        = Double
    hi Double -> Double -> Double
    forall a. Num a => a -> a -> a
    - Double
    lo
    {-# INLINABLE  kde #-}
    {-# SPECIAlIZE kde :: Int -> U.Vector Double -> (U.Vector Double, U.Vector Double) #-}
    {-# SPECIAlIZE kde :: Int -> V.Vector Double -> (V.Vector Double, V.Vector Double) #-}

    -- | Gaussian kernel density estimator for one-dimensional data, using
    -- the method of Botev et al.
    --
    -- The result is a pair of vectors, containing:
    --
    -- * The coordinates of each mesh point.
    --
    -- * Density estimates at each mesh point.
    kde_ :: (G.Vector v CD, G.Vector v Double, G.Vector v Int)
         => Int
         -- ^ The number of mesh points to use in the uniform discretization
         -- of the interval @(min,max)@.  If this value is not a power of
         -- two, then it is rounded up to the next power of two.
         -> Double
         -- ^ Lower bound (@min@) of the mesh range.
         -> Double
         -- ^ Upper bound (@max@) of the mesh range.
         -> v Double
         -> (v Double, v Double)
    kde_ :: forall (v :: * -> *).
    (Vector v CD, Vector v Double, Vector v Int) =>
    Int -> Double -> Double -> v Double -> (v Double, v Double)
    kde_ Int
    n0 Double
    min Double
    max v Double
    xs
      | v Double -> Bool
    forall (v :: * -> *) a. Vector v a => v a -> Bool
    G.null v Double
    xs = [Char] -> (v Double, v Double)
    forall a. HasCallStack => [Char] -> a
    error [Char]
    "Statistics.KernelDensity.kde: empty sample"
      | Int
    n0 Int -> Int -> Bool
    forall a. Ord a => a -> a -> Bool
    <= Int
    1   = [Char] -> (v Double, v Double)
    forall a. HasCallStack => [Char] -> a
    error [Char]
    "Statistics.KernelDensity.kde: invalid number of points"
      | Bool
    otherwise = (v Double
    mesh, v Double
    density)
      where
        mesh :: v Double
    mesh = Int -> (Int -> Double) -> v Double
    forall (v :: * -> *) a. Vector v a => Int -> (Int -> a) -> v a
    G.generate Int
    ni ((Int -> Double) -> v Double) -> (Int -> Double) -> v Double
    forall a b. (a -> b) -> a -> b
    $ \Int
    z -> Double
    min Double -> Double -> Double
    forall a. Num a => a -> a -> a
    + (Double
    d Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Int -> Double
    forall a b. (Integral a, Num b) => a -> b
    fromIntegral Int
    z)
            where d :: Double
    d = Double
    r Double -> Double -> Double
    forall a. Fractional a => a -> a -> a
    / (Double
    nDouble -> Double -> Double
    forall a. Num a => a -> a -> a
    -Double
    1)
        density :: v Double
    density = (Double -> Double) -> v Double -> v Double
    forall (v :: * -> *) a b.
    (Vector v a, Vector v b) =>
    (a -> b) -> v a -> v b
    G.map (Double -> Double -> Double
    forall a. Fractional a => a -> a -> a
    /(Double
    2 Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double
    r)) (v Double -> v Double)
    -> (v Double -> v Double) -> v Double -> v Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . v Double -> v Double
    forall (v :: * -> *).
    (Vector v CD, Vector v Double) =>
    v Double -> v Double
    idct (v Double -> v Double) -> v Double -> v Double
    forall a b. (a -> b) -> a -> b
    $ (Double -> Double -> Double) -> v Double -> v Double -> v Double
    forall (v :: * -> *) a b c.
    (Vector v a, Vector v b, Vector v c) =>
    (a -> b -> c) -> v a -> v b -> v c
    G.zipWith Double -> Double -> Double
    f v Double
    a (Double -> Double -> v Double
    forall (v :: * -> *) a. (Vector v a, Enum a) => a -> a -> v a
    G.enumFromTo Double
    0 (Double
    nDouble -> Double -> Double
    forall a. Num a => a -> a -> a
    -Double
    1))
          where f :: Double -> Double -> Double
    f Double
    b Double
    z = Double
    b Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double -> Double
    forall a. Floating a => a -> a
    exp (Double -> Double
    forall {a}. Num a => a -> a
    sqr Double
    z Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double -> Double
    forall {a}. Num a => a -> a
    sqr Double
    forall a. Floating a => a
    pi Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double
    t_star Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * (-Double
    0.5))
        !n :: Double
    n  = Int -> Double
    forall a b. (Integral a, Num b) => a -> b
    fromIntegral Int
    ni
        !ni :: Int
    ni = Int -> Int
    nextHighestPowerOfTwo Int
    n0
        !r :: Double
    r  = Double
    max Double -> Double -> Double
    forall a. Num a => a -> a -> a
    - Double
    min
        a :: v Double
    a   = v Double -> v Double
    forall (v :: * -> *).
    (Vector v CD, Vector v Double, Vector v Int) =>
    v Double -> v Double
    dct (v Double -> v Double)
    -> (v Double -> v Double) -> v Double -> v Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . (Double -> Double) -> v Double -> v Double
    forall (v :: * -> *) a b.
    (Vector v a, Vector v b) =>
    (a -> b) -> v a -> v b
    G.map (Double -> Double -> Double
    forall a. Fractional a => a -> a -> a
    / v Double -> Double
    forall (v :: * -> *). Vector v Double => v Double -> Double
    sum v Double
    h) (v Double -> v Double) -> v Double -> v Double
    forall a b. (a -> b) -> a -> b
    $ v Double
    h
            where h :: v Double
    h = (Double -> Double) -> v Double -> v Double
    forall (v :: * -> *) a b.
    (Vector v a, Vector v b) =>
    (a -> b) -> v a -> v b
    G.map (Double -> Double -> Double
    forall a. Fractional a => a -> a -> a
    / Double
    len) (v Double -> v Double) -> v Double -> v Double
    forall a b. (a -> b) -> a -> b
    $ Int -> Double -> Double -> v Double -> v Double
    forall b a (v0 :: * -> *) (v1 :: * -> *).
    (Num b, RealFrac a, Vector v0 a, Vector v1 b) =>
    Int -> a -> a -> v0 a -> v1 b
    histogram_ Int
    ni Double
    min Double
    max v Double
    xs
        !len :: Double
    len    = Int -> Double
    forall a b. (Integral a, Num b) => a -> b
    fromIntegral (v Double -> Int
    forall (v :: * -> *) a. Vector v a => v a -> Int
    G.length v Double
    xs)
        !t_star :: Double
    t_star = Double -> Root Double -> Double
    forall a. a -> Root a -> a
    fromRoot (Double
    0.28 Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double
    len Double -> Double -> Double
    forall a. Floating a => a -> a -> a
    ** (-Double
    0.4)) (Root Double -> Double)
    -> ((Double -> Double) -> Root Double)
    -> (Double -> Double)
    -> Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . RiddersParam
    -> (Double, Double) -> (Double -> Double) -> Root Double
    ridders RiddersParam
    forall a. Default a => a
    def{ riddersTol :: Tolerance
    riddersTol = Double -> Tolerance
    AbsTol Double
    1e-14 } (Double
    0,Double
    0.1)
                ((Double -> Double) -> Double) -> (Double -> Double) -> Double
    forall a b. (a -> b) -> a -> b
    $ \Double
    x -> Double
    x Double -> Double -> Double
    forall a. Num a => a -> a -> a
    - (Double
    len Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * (Double
    2 Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double -> Double
    forall a. Floating a => a -> a
    sqrt Double
    forall a. Floating a => a
    pi) Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double -> Double -> Double
    go Double
    6 (Double -> Double -> Double
    f Double
    7 Double
    x)) Double -> Double -> Double
    forall a. Floating a => a -> a -> a
    ** (-Double
    0.4)
          where
            f :: Double -> Double -> Double
    f Double
    q Double
    t = Double
    2 Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double
    forall a. Floating a => a
    pi Double -> Double -> Double
    forall a. Floating a => a -> a -> a
    ** (Double
    qDouble -> Double -> Double
    forall a. Num a => a -> a -> a
    *Double
    2) Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * v Double -> Double
    forall (v :: * -> *). Vector v Double => v Double -> Double
    sum ((Double -> Double -> Double) -> v Double -> v Double -> v Double
    forall (v :: * -> *) a b c.
    (Vector v a, Vector v b, Vector v c) =>
    (a -> b -> c) -> v a -> v b -> v c
    G.zipWith Double -> Double -> Double
    g v Double
    iv v Double
    a2v)
              where g :: Double -> Double -> Double
    g Double
    i Double
    a2 = Double
    i Double -> Double -> Double
    forall a. Floating a => a -> a -> a
    ** Double
    q Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double
    a2 Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double -> Double
    forall a. Floating a => a -> a
    exp ((-Double
    i) Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double -> Double
    forall {a}. Num a => a -> a
    sqr Double
    forall a. Floating a => a
    pi Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double
    t)
                    a2v :: v Double
    a2v = (Double -> Double) -> v Double -> v Double
    forall (v :: * -> *) a b.
    (Vector v a, Vector v b) =>
    (a -> b) -> v a -> v b
    G.map (Double -> Double
    forall {a}. Num a => a -> a
    sqr (Double -> Double) -> (Double -> Double) -> Double -> Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . (Double -> Double -> Double
    forall a. Num a => a -> a -> a
    *Double
    0.5)) (v Double -> v Double) -> v Double -> v Double
    forall a b. (a -> b) -> a -> b
    $ v Double -> v Double
    forall (v :: * -> *) a. Vector v a => v a -> v a
    G.tail v Double
    a
                    iv :: v Double
    iv = (Double -> Double) -> v Double -> v Double
    forall (v :: * -> *) a b.
    (Vector v a, Vector v b) =>
    (a -> b) -> v a -> v b
    G.map Double -> Double
    forall {a}. Num a => a -> a
    sqr (v Double -> v Double) -> v Double -> v Double
    forall a b. (a -> b) -> a -> b
    $ Double -> Double -> v Double
    forall (v :: * -> *) a. (Vector v a, Enum a) => a -> a -> v a
    G.enumFromTo Double
    1 (Double
    nDouble -> Double -> Double
    forall a. Num a => a -> a -> a
    -Double
    1)
            go :: Double -> Double -> Double
    go Double
    s !Double
    h | Double
    s Double -> Double -> Bool
    forall a. Eq a => a -> a -> Bool
    == Double
    1    = Double
    h
                    | Bool
    otherwise = Double -> Double -> Double
    go (Double
    sDouble -> Double -> Double
    forall a. Num a => a -> a -> a
    -Double
    1) (Double -> Double -> Double
    f Double
    s Double
    time)
              where time :: Double
    time  = (Double
    2 Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double
    const Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double
    k0 Double -> Double -> Double
    forall a. Fractional a => a -> a -> a
    / Double
    len Double -> Double -> Double
    forall a. Fractional a => a -> a -> a
    / Double
    h) Double -> Double -> Double
    forall a. Floating a => a -> a -> a
    ** (Double
    2 Double -> Double -> Double
    forall a. Fractional a => a -> a -> a
    / (Double
    3 Double -> Double -> Double
    forall a. Num a => a -> a -> a
    + Double
    2 Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double
    s))
                    const :: Double
    const = (Double
    1 Double -> Double -> Double
    forall a. Num a => a -> a -> a
    + Double
    0.5 Double -> Double -> Double
    forall a. Floating a => a -> a -> a
    ** (Double
    sDouble -> Double -> Double
    forall a. Num a => a -> a -> a
    +Double
    0.5)) Double -> Double -> Double
    forall a. Fractional a => a -> a -> a
    / Double
    3
                    k0 :: Double
    k0    = Vector Double -> Double
    forall a. (Unbox a, Num a) => Vector a -> a
    U.product (Double -> Double -> Double -> Vector Double
    forall (v :: * -> *) a. (Vector v a, Enum a) => a -> a -> a -> v a
    G.enumFromThenTo Double
    1 Double
    3 (Double
    2Double -> Double -> Double
    forall a. Num a => a -> a -> a
    *Double
    sDouble -> Double -> Double
    forall a. Num a => a -> a -> a
    -Double
    1)) Double -> Double -> Double
    forall a. Fractional a => a -> a -> a
    / Double
    m_sqrt_2_pi
        sqr :: a -> a
    sqr a
    x = a
    x a -> a -> a
    forall a. Num a => a -> a -> a
    * a
    x
    {-# INLINABLE  kde_ #-}
    {-# SPECIAlIZE kde_ :: Int -> Double -> Double -> U.Vector Double -> (U.Vector Double, U.Vector Double) #-}
    {-# SPECIAlIZE kde_ :: Int -> Double -> Double -> V.Vector Double -> (V.Vector Double, V.Vector Double) #-}

    -- $references
    --
    -- Botev. Z.I., Grotowski J.F., Kroese D.P. (2010). Kernel density
    -- estimation via diffusion. /Annals of Statistics/
    -- 38(5):2916&#8211;2957. <http://arxiv.org/pdf/1011.2602>