{-# LANGUAGE DeriveDataTypeable, DeriveGeneric, FlexibleContexts #-}
module Statistics.Sample.KernelDensity.Simple
{-# DEPRECATED "Use Statistics.Sample.KernelDensity instead." #-}
(
epanechnikovPDF
, gaussianPDF
, Points(..)
, choosePoints
, Bandwidth
, bandwidth
, epanechnikovBW
, gaussianBW
, Kernel
, epanechnikovKernel
, gaussianKernel
, estimatePDF
, simplePDF
) where
import Data.Aeson (FromJSON, ToJSON)
import Data.Binary (Binary(..))
import Data.Data (Data, Typeable)
import Data.Vector.Binary ()
import GHC.Generics (Generic)
import Numeric.MathFunctions.Constants (m_1_sqrt_2, m_2_sqrt_pi)
import Prelude hiding (sum)
import Statistics.Function (minMax)
import Statistics.Sample (stdDev)
import Statistics.Sample.Internal (sum)
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Unboxed as U
newtype Points = Points {
Points -> Vector Double
fromPoints :: U.Vector Double
} deriving (Points -> Points -> Bool
(Points -> Points -> Bool)
-> (Points -> Points -> Bool) -> Eq Points
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: Points -> Points -> Bool
$c/= :: Points -> Points -> Bool
== :: Points -> Points -> Bool
$c== :: Points -> Points -> Bool
Eq, ReadPrec [Points]
ReadPrec Points
Int -> ReadS Points
ReadS [Points]
(Int -> ReadS Points)
-> ReadS [Points]
-> ReadPrec Points
-> ReadPrec [Points]
-> Read Points
forall a.
(Int -> ReadS a)
-> ReadS [a] -> ReadPrec a -> ReadPrec [a] -> Read a
readListPrec :: ReadPrec [Points]
$creadListPrec :: ReadPrec [Points]
readPrec :: ReadPrec Points
$creadPrec :: ReadPrec Points
readList :: ReadS [Points]
$creadList :: ReadS [Points]
readsPrec :: Int -> ReadS Points
$creadsPrec :: Int -> ReadS Points
Read, Int -> Points -> ShowS
[Points] -> ShowS
Points -> String
(Int -> Points -> ShowS)
-> (Points -> String) -> ([Points] -> ShowS) -> Show Points
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [Points] -> ShowS
$cshowList :: [Points] -> ShowS
show :: Points -> String
$cshow :: Points -> String
showsPrec :: Int -> Points -> ShowS
$cshowsPrec :: Int -> Points -> ShowS
Show, Typeable, Typeable Points
DataType
Constr
Typeable Points
-> (forall (c :: * -> *).
(forall d b. Data d => c (d -> b) -> d -> c b)
-> (forall g. g -> c g) -> Points -> c Points)
-> (forall (c :: * -> *).
(forall b r. Data b => c (b -> r) -> c r)
-> (forall r. r -> c r) -> Constr -> c Points)
-> (Points -> Constr)
-> (Points -> DataType)
-> (forall (t :: * -> *) (c :: * -> *).
Typeable t =>
(forall d. Data d => c (t d)) -> Maybe (c Points))
-> (forall (t :: * -> * -> *) (c :: * -> *).
Typeable t =>
(forall d e. (Data d, Data e) => c (t d e)) -> Maybe (c Points))
-> ((forall b. Data b => b -> b) -> Points -> Points)
-> (forall r r'.
(r -> r' -> r)
-> r -> (forall d. Data d => d -> r') -> Points -> r)
-> (forall r r'.
(r' -> r -> r)
-> r -> (forall d. Data d => d -> r') -> Points -> r)
-> (forall u. (forall d. Data d => d -> u) -> Points -> [u])
-> (forall u. Int -> (forall d. Data d => d -> u) -> Points -> u)
-> (forall (m :: * -> *).
Monad m =>
(forall d. Data d => d -> m d) -> Points -> m Points)
-> (forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> Points -> m Points)
-> (forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> Points -> m Points)
-> Data Points
Points -> DataType
Points -> Constr
(forall b. Data b => b -> b) -> Points -> Points
(forall d b. Data d => c (d -> b) -> d -> c b)
-> (forall g. g -> c g) -> Points -> c Points
(forall b r. Data b => c (b -> r) -> c r)
-> (forall r. r -> c r) -> Constr -> c Points
forall a.
Typeable a
-> (forall (c :: * -> *).
(forall d b. Data d => c (d -> b) -> d -> c b)
-> (forall g. g -> c g) -> a -> c a)
-> (forall (c :: * -> *).
(forall b r. Data b => c (b -> r) -> c r)
-> (forall r. r -> c r) -> Constr -> c a)
-> (a -> Constr)
-> (a -> DataType)
-> (forall (t :: * -> *) (c :: * -> *).
Typeable t =>
(forall d. Data d => c (t d)) -> Maybe (c a))
-> (forall (t :: * -> * -> *) (c :: * -> *).
Typeable t =>
(forall d e. (Data d, Data e) => c (t d e)) -> Maybe (c a))
-> ((forall b. Data b => b -> b) -> a -> a)
-> (forall r r'.
(r -> r' -> r) -> r -> (forall d. Data d => d -> r') -> a -> r)
-> (forall r r'.
(r' -> r -> r) -> r -> (forall d. Data d => d -> r') -> a -> r)
-> (forall u. (forall d. Data d => d -> u) -> a -> [u])
-> (forall u. Int -> (forall d. Data d => d -> u) -> a -> u)
-> (forall (m :: * -> *).
Monad m =>
(forall d. Data d => d -> m d) -> a -> m a)
-> (forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> a -> m a)
-> (forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> a -> m a)
-> Data a
forall u. Int -> (forall d. Data d => d -> u) -> Points -> u
forall u. (forall d. Data d => d -> u) -> Points -> [u]
forall r r'.
(r -> r' -> r) -> r -> (forall d. Data d => d -> r') -> Points -> r
forall r r'.
(r' -> r -> r) -> r -> (forall d. Data d => d -> r') -> Points -> r
forall (m :: * -> *).
Monad m =>
(forall d. Data d => d -> m d) -> Points -> m Points
forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> Points -> m Points
forall (c :: * -> *).
(forall b r. Data b => c (b -> r) -> c r)
-> (forall r. r -> c r) -> Constr -> c Points
forall (c :: * -> *).
(forall d b. Data d => c (d -> b) -> d -> c b)
-> (forall g. g -> c g) -> Points -> c Points
forall (t :: * -> *) (c :: * -> *).
Typeable t =>
(forall d. Data d => c (t d)) -> Maybe (c Points)
forall (t :: * -> * -> *) (c :: * -> *).
Typeable t =>
(forall d e. (Data d, Data e) => c (t d e)) -> Maybe (c Points)
$cPoints :: Constr
$tPoints :: DataType
gmapMo :: (forall d. Data d => d -> m d) -> Points -> m Points
$cgmapMo :: forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> Points -> m Points
gmapMp :: (forall d. Data d => d -> m d) -> Points -> m Points
$cgmapMp :: forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> Points -> m Points
gmapM :: (forall d. Data d => d -> m d) -> Points -> m Points
$cgmapM :: forall (m :: * -> *).
Monad m =>
(forall d. Data d => d -> m d) -> Points -> m Points
gmapQi :: Int -> (forall d. Data d => d -> u) -> Points -> u
$cgmapQi :: forall u. Int -> (forall d. Data d => d -> u) -> Points -> u
gmapQ :: (forall d. Data d => d -> u) -> Points -> [u]
$cgmapQ :: forall u. (forall d. Data d => d -> u) -> Points -> [u]
gmapQr :: (r' -> r -> r) -> r -> (forall d. Data d => d -> r') -> Points -> r
$cgmapQr :: forall r r'.
(r' -> r -> r) -> r -> (forall d. Data d => d -> r') -> Points -> r
gmapQl :: (r -> r' -> r) -> r -> (forall d. Data d => d -> r') -> Points -> r
$cgmapQl :: forall r r'.
(r -> r' -> r) -> r -> (forall d. Data d => d -> r') -> Points -> r
gmapT :: (forall b. Data b => b -> b) -> Points -> Points
$cgmapT :: (forall b. Data b => b -> b) -> Points -> Points
dataCast2 :: (forall d e. (Data d, Data e) => c (t d e)) -> Maybe (c Points)
$cdataCast2 :: forall (t :: * -> * -> *) (c :: * -> *).
Typeable t =>
(forall d e. (Data d, Data e) => c (t d e)) -> Maybe (c Points)
dataCast1 :: (forall d. Data d => c (t d)) -> Maybe (c Points)
$cdataCast1 :: forall (t :: * -> *) (c :: * -> *).
Typeable t =>
(forall d. Data d => c (t d)) -> Maybe (c Points)
dataTypeOf :: Points -> DataType
$cdataTypeOf :: Points -> DataType
toConstr :: Points -> Constr
$ctoConstr :: Points -> Constr
gunfold :: (forall b r. Data b => c (b -> r) -> c r)
-> (forall r. r -> c r) -> Constr -> c Points
$cgunfold :: forall (c :: * -> *).
(forall b r. Data b => c (b -> r) -> c r)
-> (forall r. r -> c r) -> Constr -> c Points
gfoldl :: (forall d b. Data d => c (d -> b) -> d -> c b)
-> (forall g. g -> c g) -> Points -> c Points
$cgfoldl :: forall (c :: * -> *).
(forall d b. Data d => c (d -> b) -> d -> c b)
-> (forall g. g -> c g) -> Points -> c Points
$cp1Data :: Typeable Points
Data, (forall x. Points -> Rep Points x)
-> (forall x. Rep Points x -> Points) -> Generic Points
forall x. Rep Points x -> Points
forall x. Points -> Rep Points x
forall a.
(forall x. a -> Rep a x) -> (forall x. Rep a x -> a) -> Generic a
$cto :: forall x. Rep Points x -> Points
$cfrom :: forall x. Points -> Rep Points x
Generic)
instance FromJSON Points
instance ToJSON Points
instance Binary Points where
get :: Get Points
get = (Vector Double -> Points) -> Get (Vector Double) -> Get Points
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap Vector Double -> Points
Points Get (Vector Double)
forall t. Binary t => Get t
get
put :: Points -> Put
put = Vector Double -> Put
forall t. Binary t => t -> Put
put (Vector Double -> Put)
-> (Points -> Vector Double) -> Points -> Put
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Points -> Vector Double
fromPoints
epanechnikovBW :: Double -> Bandwidth
epanechnikovBW :: Double -> Double
epanechnikovBW Double
n = (Double
80 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
m_2_sqrt_pi)) Double -> Double -> Double
forall a. Floating a => a -> a -> a
** Double
0.2
gaussianBW :: Double -> Bandwidth
gaussianBW :: Double -> Double
gaussianBW Double
n = (Double
4 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
3)) Double -> Double -> Double
forall a. Floating a => a -> a -> a
** Double
0.2
type Bandwidth = Double
bandwidth :: G.Vector v Double =>
(Double -> Bandwidth)
-> v Double
-> Bandwidth
bandwidth :: (Double -> Double) -> v Double -> Double
bandwidth Double -> Double
kern v Double
values = v Double -> Double
forall (v :: * -> *). Vector v Double => v Double -> Double
stdDev v Double
values Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
kern (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Double) -> Int -> Double
forall a b. (a -> b) -> a -> b
$ v Double -> Int
forall (v :: * -> *) a. Vector v a => v a -> Int
G.length v Double
values)
choosePoints :: G.Vector v Double =>
Int
-> Double
-> v Double
-> Points
choosePoints :: Int -> Double -> v Double -> Points
choosePoints Int
n Double
h v Double
sample = Vector Double -> Points
Points (Vector Double -> Points)
-> (Vector Int -> Vector Double) -> Vector Int -> Points
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int -> Double) -> Vector Int -> Vector Double
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
U.map Int -> Double
forall a. Integral a => a -> Double
f (Vector Int -> Points) -> Vector Int -> Points
forall a b. (a -> b) -> a -> b
$ Int -> Int -> Vector Int
forall a. (Unbox a, Enum a) => a -> a -> Vector a
U.enumFromTo Int
0 Int
n'
where lo :: Double
lo = Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
h
hi :: Double
hi = Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
h
(Double
a, Double
z) = v Double -> (Double, Double)
forall (v :: * -> *).
Vector v Double =>
v Double -> (Double, Double)
minMax v Double
sample
d :: Double
d = (Double
hi Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
lo) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n'
f :: a -> Double
f a
i = Double
lo Double -> Double -> Double
forall a. Num a => a -> a -> a
+ a -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
i Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d
n' :: Int
n' = Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1
type Kernel = Double
-> Double
-> Double
-> Double
-> Double
epanechnikovKernel :: Kernel
epanechnikovKernel :: Kernel
epanechnikovKernel Double
f Double
h Double
p Double
v
| Double -> Double
forall a. Num a => a -> a
abs Double
u Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
1 = Double
f Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u)
| Bool
otherwise = Double
0
where u :: Double
u = (Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
p) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
h Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
0.75)
gaussianKernel :: Kernel
gaussianKernel :: Kernel
gaussianKernel Double
f Double
h Double
p Double
v = Double -> Double
forall a. Floating a => a -> a
exp (-Double
0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
g
where u :: Double
u = (Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
p) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
h
g :: Double
g = Double
f Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
m_2_sqrt_pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
m_1_sqrt_2
estimatePDF :: G.Vector v Double =>
Kernel
-> Bandwidth
-> v Double
-> Points
-> U.Vector Double
estimatePDF :: Kernel -> Double -> v Double -> Points -> Vector Double
estimatePDF Kernel
kernel Double
h v Double
sample
| Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
2 = String -> Points -> Vector Double
forall a. String -> a
errorShort String
"estimatePDF"
| Bool
otherwise = (Double -> Double) -> Vector Double -> Vector Double
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
U.map Double -> Double
k (Vector Double -> Vector Double)
-> (Points -> Vector Double) -> Points -> Vector Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Points -> Vector Double
fromPoints
where
k :: Double -> Double
k Double
p = v Double -> Double
forall (v :: * -> *). Vector v Double => v Double -> Double
sum (v Double -> Double)
-> (v Double -> v Double) -> v Double -> 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 (Kernel
kernel Double
f Double
h Double
p) (v Double -> Double) -> v Double -> Double
forall a b. (a -> b) -> a -> b
$ v Double
sample
f :: Double
f = Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
h Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n)
n :: Int
n = v Double -> Int
forall (v :: * -> *) a. Vector v a => v a -> Int
G.length v Double
sample
{-# INLINE estimatePDF #-}
simplePDF :: G.Vector v Double =>
(Double -> Double)
-> Kernel
-> Double
-> Int
-> v Double
-> (Points, U.Vector Double)
simplePDF :: (Double -> Double)
-> Kernel -> Double -> Int -> v Double -> (Points, Vector Double)
simplePDF Double -> Double
fbw Kernel
fpdf Double
k Int
numPoints v Double
sample =
(Points
points, Kernel -> Double -> v Double -> Points -> Vector Double
forall (v :: * -> *).
Vector v Double =>
Kernel -> Double -> v Double -> Points -> Vector Double
estimatePDF Kernel
fpdf Double
bw v Double
sample Points
points)
where points :: Points
points = Int -> Double -> v Double -> Points
forall (v :: * -> *).
Vector v Double =>
Int -> Double -> v Double -> Points
choosePoints Int
numPoints (Double
bwDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
k) v Double
sample
bw :: Double
bw = (Double -> Double) -> v Double -> Double
forall (v :: * -> *).
Vector v Double =>
(Double -> Double) -> v Double -> Double
bandwidth Double -> Double
fbw v Double
sample
{-# INLINE simplePDF #-}
epanechnikovPDF :: G.Vector v Double =>
Int
-> v Double
-> (Points, U.Vector Double)
epanechnikovPDF :: Int -> v Double -> (Points, Vector Double)
epanechnikovPDF = (Double -> Double)
-> Kernel -> Double -> Int -> v Double -> (Points, Vector Double)
forall (v :: * -> *).
Vector v Double =>
(Double -> Double)
-> Kernel -> Double -> Int -> v Double -> (Points, Vector Double)
simplePDF Double -> Double
epanechnikovBW Kernel
epanechnikovKernel Double
1
gaussianPDF :: G.Vector v Double =>
Int
-> v Double
-> (Points, U.Vector Double)
gaussianPDF :: Int -> v Double -> (Points, Vector Double)
gaussianPDF = (Double -> Double)
-> Kernel -> Double -> Int -> v Double -> (Points, Vector Double)
forall (v :: * -> *).
Vector v Double =>
(Double -> Double)
-> Kernel -> Double -> Int -> v Double -> (Points, Vector Double)
simplePDF Double -> Double
gaussianBW Kernel
gaussianKernel Double
3
errorShort :: String -> a
errorShort :: String -> a
errorShort String
func = String -> a
forall a. HasCallStack => String -> a
error (String
"Statistics.KernelDensity." String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
func String -> ShowS
forall a. [a] -> [a] -> [a]
++
String
": at least two points required")