{-# LANGUAGE BangPatterns, FlexibleContexts, UnboxedTuples #-}
module Statistics.Sample.KernelDensity
(
kde
, kde_
) 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
kde :: (G.Vector v CD, G.Vector v Double, G.Vector v Int)
=> Int
-> v Double -> (v Double, v Double)
kde :: Int -> v Double -> (v Double, v Double)
kde n0 :: Int
n0 xs :: 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
/ 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
/ 10) v Double
xs
where
(lo :: Double
lo,hi :: 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
<= 1 = 1
| Double
lo Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
hi = 1
| 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) #-}
kde_ :: (G.Vector v CD, G.Vector v Double, G.Vector v Int)
=> Int
-> Double
-> Double
-> v Double
-> (v Double, v Double)
kde_ :: Int -> Double -> Double -> v Double -> (v Double, v Double)
kde_ n0 :: Int
n0 min :: Double
min max :: Double
max xs :: 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 "Statistics.KernelDensity.kde: empty sample"
| Int
n0 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= 1 = [Char] -> (v Double, v Double)
forall a. HasCallStack => [Char] -> a
error "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
$ \z :: 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
-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
/(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 0 (Double
nDouble -> Double -> Double
forall a. Num a => a -> a -> a
-1))
where f :: Double -> Double -> Double
f b :: Double
b z :: 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
* (-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 (0.28 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
len Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (-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 1e-14 } (0,0.1)
((Double -> Double) -> Double) -> (Double -> Double) -> Double
forall a b. (a -> b) -> a -> b
$ \x :: 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
* (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 6 (Double -> Double -> Double
f 7 Double
x)) Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (-0.4)
where
f :: Double -> Double -> Double
f q :: Double
q t :: Double
t = 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
*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 i :: Double
i a2 :: 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
*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 1 (Double
nDouble -> Double -> Double
forall a. Num a => a -> a -> a
-1)
go :: Double -> Double -> Double
go s :: Double
s !Double
h | Double
s Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== 1 = Double
h
| Bool
otherwise = Double -> Double -> Double
go (Double
sDouble -> Double -> Double
forall a. Num a => a -> a -> a
-1) (Double -> Double -> Double
f Double
s Double
time)
where time :: Double
time = (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
** (2 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (3 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
s))
const :: Double
const = (1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 0.5 Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (Double
sDouble -> Double -> Double
forall a. Num a => a -> a -> a
+0.5)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ 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 1 3 (2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
sDouble -> Double -> Double
forall a. Num a => a -> a -> a
-1)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
m_sqrt_2_pi
sqr :: a -> a
sqr x :: 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) #-}