HLearn-distributions 0.0.1.3 → 0.1.0.0

raw patch · 10 files changed

+573/−495 lines, 10 filesdep +deepseqdep +processdep −math-functionsdep −normaldistributiondep ~ConstraintKindsdep ~HLearn-algebradep ~logfloatPVP ok

version bump matches the API change (PVP)

Dependencies added: deepseq, process

Dependencies removed: math-functions, normaldistribution

Dependency ranges changed: ConstraintKinds, HLearn-algebra, logfloat

API changes (from Hackage documentation)

Files

• HLearn-distributions.cabal +12/−8
• src/HLearn/Gnuplot/Distributions.hs +79/−0
• src/HLearn/Models/Distributions.hs +4/−2
• src/HLearn/Models/Distributions/Categorical.hs +33/−10
• src/HLearn/Models/Distributions/Common.hs +4/−7
• src/HLearn/Models/Distributions/Gaussian.lhs +0/−464
• src/HLearn/Models/Distributions/KernelDensityEstimator.hs +257/−0
• src/HLearn/Models/Distributions/Moments.hs +178/−0
• src/examples/Criterion.hs +3/−2
• src/examples/SpaceTests.hs +3/−2

@@ -1,5 +1,5 @@ Name:                HLearn-distributions-Version:             0.0.1.3+Version:             0.1.0.0 Synopsis:            Distributions for use with the HLearn library Description:         This module is used to estimate statistical distributions from data.  The focus is a clean interface inspired by algebra. Category:            Data Mining, Machine Learning, Statistics@@ -54,9 +54,11 @@  Library     Build-Depends:      -        --HLearn-algebra          >= 0.0.1,-        HLearn-algebra          ,+        HLearn-algebra          >= 0.1.0,+        ConstraintKinds         >= 0.0.1,         base                    >= 3 && < 5,+        +        deepseq                 >= 1.3.0.1,         list-extras             >= 0.4.1,         containers              >= 0.5,         statistics              >= 0.10.2,@@ -66,10 +68,8 @@         vector-th-unbox         >= 0.2,                  -- are these really necessary?-        MonadRandom             >= 0.1.6,-        logfloat                >= 0.12.0,-        math-functions          >= 0.1.1,-        normaldistribution      >= 1.1.0+        process                 >= 1.1.0.2,+        MonadRandom             >= 0.1.6               hs-source-dirs:     src@@ -81,10 +81,14 @@         -O2          -- -fllvm     Exposed-modules:+        HLearn.Gnuplot.Distributions         HLearn.Models.Distributions         HLearn.Models.Distributions.Common         HLearn.Models.Distributions.Categorical-        HLearn.Models.Distributions.Gaussian+        HLearn.Models.Distributions.KernelDensityEstimator+        HLearn.Models.Distributions.Moments+        --HLearn.Models.Distributions.Normal+        --HLearn.Models.Distributions.Gaussian         --HLearn.Models.Distributions.GaussianOld         --HLearn.Models.Distributions.GaussianOld2         --HLearn.Models.Distributions.Poisson

@@ -0,0 +1,79 @@+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE TypeSynonymInstances #-}++module HLearn.Gnuplot.Distributions+    where++import HLearn.Algebra+import HLearn.Models.Distributions++import qualified Data.Vector.Unboxed as VU++import Control.Monad+import System.IO+import System.Process++data PlotParams = PlotParams+    { dataFile :: FilePath+    , gnuFile :: FilePath+    , picFile :: FilePath+    , plotPoints :: [Double]+    }++class PlottableDistribution dist where+    genPlotPoints :: dist -> [Double]++instance PlottableDistribution (KDE Double) where+    genPlotPoints (SGJust kde) = VU.toList $ samplePoints $ params kde++genPlotParams :: (PlottableDistribution dist) => String -> dist -> PlotParams+genPlotParams str dist = PlotParams+    { dataFile = str++".dat"+    , gnuFile  = str++".gnu"+    , picFile  = str++".ps"+    , plotPoints = genPlotPoints dist+    }++defPlotParams = PlotParams+    { dataFile = "hlearn-distributions.dat"+    , gnuFile = "hlearn-distributions.gnu"+    , picFile = "hlearn-distributions.ps"+    , plotPoints = [-5..5]+    }+    +plotDistribution :: (Distribution dist Double Double) => PlotParams -> dist -> IO ()+plotDistribution params dist = do+    +    -- Create data file+    putStrLn "Creating data file..."+    datah <- openFile (dataFile params) WriteMode+    forM_ (plotPoints params) $ \x -> do+--     forM_ (map (/10) [-50..50]) $ \x -> do+        hPutStrLn datah $ show (x::Double) ++ " " ++ show (pdf dist x::Double)+    hClose datah+    +    -- Create gnuplot file+    putStrLn "Plotting data"+    gnuh <- openFile (gnuFile params) WriteMode+    hPutStrLn gnuh $ "set terminal postscript \"Times-Roman\" 25"+    hPutStrLn gnuh $ "set output \"" ++ (picFile params) ++ "\""+    hPutStrLn gnuh $ "unset xtics; unset ytics; unset key"+    hPutStrLn gnuh $ "set border 0; set xzeroaxis lt 1; set yzeroaxis lt 1"+    hPutStrLn gnuh $ "plot '"++(dataFile params)++"' using 1:2 lt 1 lw 4 lc rgb '#ccccff' with filledcurves, \\"+    hPutStrLn gnuh $ "     '"++(dataFile params)++"' using 1:2 lt 1 lw 4 lc rgb '#0000ff' with lines"+    hClose gnuh+    +    -- Run gnuplot, creating picture+    system $ "gnuplot "++(gnuFile params)+    putStrLn "done."+    return ()+    +-- testL = [1,1.1,1.2,0,5,-3,2,1,2.5,2.9::Double]+-- kde = train' (KDEParams 1 (VU.fromList $ map (/10) [-50..50::Double]) (KernelBox Gaussian)) testL :: KDE Double+-- +-- pdfL = sequence_ $ do+--     x <- map (/10) [-50..50]+--     let y = pdf kde x :: Double+--     seq y (return ())+--     return $ putStrLn $ (show (x::Double)) ++ " " ++ show y

@@ -6,13 +6,15 @@ module HLearn.Models.Distributions     ( module HLearn.Models.Distributions.Common     , module HLearn.Models.Distributions.Categorical-    , module HLearn.Models.Distributions.Gaussian+    , module HLearn.Models.Distributions.KernelDensityEstimator+--     , module HLearn.Models.Distributions.Gaussian     )     where  import HLearn.Models.Distributions.Common import HLearn.Models.Distributions.Categorical-import HLearn.Models.Distributions.Gaussian+import HLearn.Models.Distributions.KernelDensityEstimator+-- import HLearn.Models.Distributions.Gaussian -- import HLearn.DataContainers -- import HLearn.Math.TypeClasses 

@@ -9,9 +9,11 @@     ( Categorical (..)     , CategoricalParams(..)     , dist2list+    , mostLikely     )     where +import Control.DeepSeq import Control.Monad.Random import Data.List import Data.List.Extras@@ -28,7 +30,7 @@  -- | The Categorical distribution takes no parameters data CategoricalParams = CategoricalParams-    deriving (Read,Show,Eq)+    deriving (Read,Show,Eq,Ord)  instance NFData CategoricalParams where     rnf x = ()@@ -36,7 +38,8 @@ instance Model CategoricalParams (Categorical label probtype) where     getparams model = CategoricalParams -instance DefaultModel CategoricalParams (Categorical label probtype) where+-- instance DefaultModel CategoricalParams (Categorical label probtype) where+instance DefaultModel CategoricalParams (Categorical Int Double) where     defparams = CategoricalParams  -------------------------------------------------------------------------------@@ -45,7 +48,7 @@ data Categorical sampletype probtype = Categorical          { pdfmap :: !(Map.Map sampletype probtype)         } -    deriving (Show,Read,Eq)+    deriving (Show,Read,Eq,Ord)  dist2list :: Categorical sampletype probtype -> [(sampletype,probtype)] dist2list (Categorical pdfmap) = Map.toList pdfmap@@ -56,9 +59,6 @@ ------------------------------------------------------------------------------- -- Training -instance (Ord label, Num probtype) => HomTrainer CategoricalParams (label,probtype) (Categorical label probtype) where-    train1dp' params (dp,w) = Categorical $ Map.singleton dp w- instance (Ord label, Num probtype) => HomTrainer CategoricalParams label (Categorical label probtype) where     train1dp' params dp = Categorical $ Map.singleton dp 1 @@ -75,7 +75,7 @@                 Just x  -> x             tot = F.foldl' (+) 0 $ pdfmap dist -    {-# INLINE cdf #-}+{-    {-# INLINE cdf #-}     cdf dist label = (Map.foldl' (+) 0 $ Map.filterWithKey (\k a -> k<=label) $ pdfmap dist)                     / (Map.foldl' (+) 0 $ pdfmap dist)                    @@ -89,7 +89,7 @@                 else go (prob-snd x) xs --     cdfInverse dist prob = argmax (cdf dist) $ Map.keys $ pdfmap dist -{-    {-# INLINE mean #-}+    {-# INLINE mean #-}     mean dist = fst $ argmax snd $ Map.toList $ pdfmap dist      {-# INLINE drawSample #-}@@ -97,10 +97,15 @@         x <- getRandomR (0,1)         return $ cdfInverse dist (x::prob) -}-        ++mostLikely :: Ord prob => Categorical label prob -> label+mostLikely dist = fst $ argmax snd $ Map.toList $ pdfmap dist++ ------------------------------------------------------------------------------- -- Algebra +instance (Ord label, Num probtype{-, NFData probtype-}) => Abelian (Categorical label probtype) instance (Ord label, Num probtype{-, NFData probtype-}) => Semigroup (Categorical label probtype) where     (<>) !d1 !d2 = {-deepseq res $-} Categorical $ res         where@@ -113,4 +118,22 @@     mempty = Categorical Map.empty     mappend = (<>) -instance (Ord label, Num probtype) => Group (Categorical label probtype)+-- instance (Ord label, Num probtype) => Group (Categorical label probtype)++instance (Ord label, Num probtype) => LeftModule probtype (Categorical label probtype)+instance (Ord label, Num probtype) => LeftOperator probtype (Categorical label probtype) where+    p .* (Categorical pdf) = Categorical $ Map.map (*p) pdf++instance (Ord label, Num probtype) => RightModule probtype (Categorical label probtype)+instance (Ord label, Num probtype) => RightOperator probtype (Categorical label probtype) where+    (*.) = flip (.*)++-------------------------------------------------------------------------------+-- Morphisms++instance +    ( Ord label+    , Num probtype+    ) => Morphism (Categorical label probtype) FreeModParams (FreeMod probtype label) +        where+    Categorical pdf $> FreeModParams = FreeMod pdf

@@ -3,22 +3,19 @@  module HLearn.Models.Distributions.Common     ( Distribution(..)-    , module Data.Number.LogFloat     )     where -import Control.Monad.Random-import Data.Number.LogFloat hiding (log) import HLearn.Algebra  ------------------------------------------------------------------------------- -- Distribution      -- | We use the same class for both discrete and continuous distributions.  Unfortunately, we cannot use the type classes from the 'statistics' package because we require more generalilty.-class Distribution dist sample prob | dist -> sample where-    pdf :: dist -> sample -> prob -    cdf :: dist -> sample -> prob -    cdfInverse :: dist -> prob -> sample+class Distribution dist dp prob | dist -> dp where+    pdf :: dist -> dp -> prob +--     cdf :: dist -> dp -> prob +--     cdfInverse :: dist -> prob -> dp  --     mean :: dist -> sample --     drawSample :: (RandomGen g) => dist -> Rand g sample

@@ -1,464 +0,0 @@-\documentclass{article}-%include polycode.fmt--\usepackage{amsmath,amsthm,amsfonts,amssymb,graphicx,color}--\newcommand{\todo}[1]{\textcolor{red}{[TODO: #1]}}-%\newcommand{\todo}[1]{[TODO: #1]}--\newtheorem{theorem}{Theorem}-\newtheorem{claim}{Claim}-\newtheorem{lemma}{Lemma}-\newtheorem{corollary}{Corollary}--\newtheoremstyle{namedtheorem}{}{}{\itshape}{}{\bfseries}{.}{.5em}{\thmnote{#3 }#1}-\theoremstyle{namedtheorem}-\newtheorem*{namedtheorem}{Theorem}-\newtheorem*{namedcorollary}{Corollary}--\newtheoremstyle{nameddefinition2}{}{}{}{}{\bfseries}{}{.5em}{#1\thmnote{#3}.}-\theoremstyle{nameddefinition2}-\newtheorem*{nnote}{}---\long\def\ignore#1{}--%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%-% Document info-% --\title{Algebraic Properties of the Gaussian Distribution}--\author{-Michael Izbicki-}--% -% Document-% --\begin{document}-\maketitle--%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%-% lhs2TeX Formating-% --%format m1 = m"_1"-%format m2 = m"_2"-%format m1a = m"_{1a}"-%format m2a = m"_{2a}"-%format na = n"_a"-%format m1b = m"_{1b}"-%format m2b = m"_{2b}"-%format nb = n"_b"-%format m1' = m"_1''"-%format m2' = m"_2''"-%format n' = n"''"--%format g_a = g"_a"-%format g_b = g"_b"-%format g_c = g"_c"--%format ^ = "^"-%format <> = "\diamond"-%format <> = "\diamond"-%format Infinity = "\infty"--%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%-% Code Header-%--\ignore{-\begin{code}-{-# LANGUAGE FlexibleInstances #-}-{-# LANGUAGE FlexibleContexts #-}-{-# LANGUAGE MultiParamTypeClasses #-}-{-# LANGUAGE UndecidableInstances #-}-{-# LANGUAGE OverlappingInstances #-}-{-# LANGUAGE BangPatterns #-}--{-# LANGUAGE TemplateHaskell #-}-{-# LANGUAGE TypeFamilies #-}---- | The Gaussian distribution is an instance of 'HomTrainer.'  For examples of how to use this type, and the math behind it, see: <http://izbicki.me/blog/gausian-distributions-are-monoids>.--module HLearn.Models.Distributions.Gaussian-    ( Gaussian (..)-    , GaussianParams (..)-    )-    where-    -import Control.Monad-import Control.Monad.Random-import Data.Random.Normal-import Debug.Trace-import Numeric.SpecFunctions (logFactorial)-import Test.QuickCheck--import Data.Vector.Generic.Base-import Data.Vector.Generic.Mutable-import qualified Data.Vector.Unboxed as U-import Data.Vector.Unboxed.Deriving--import qualified Statistics.Distribution as D-import qualified Statistics.Distribution.Normal as N--import HLearn.Algebra-import HLearn.Models.Distributions.Common--{-# INLINE fi #-}-fi :: (Integral a, Num b) => a -> b-fi=fromIntegral--\end{code}-}--%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%-% Document Begins Here-%--\section{Introduction}--\subsection{About this document}-This document covers two things.  First, it covers the theory of using algebra to manipulate Gaussian distributions.  To my knowledge, this is fairly unexplored territory.  Second, this document covers the design decisions made in HLearn's implementation of this algebra.  This document is specifically not meant to be a reference manual for using the HLearn library.  That's what the Haddock documentation is for.  -Finally, the HLearn library was written using literate Haskell, so this file is actually also the source code to the library!--\subsection{About the Gaussian distribution}-The Gaussian distribution (also called the Normal distribution) is one of the most common distributions in statistics and machine learning.--\section{Data Types}-A Gaussian distribution is uniquely determined from its mean and variance.  In order to estimate these quantities from a sample, we must also know the number of data points in the sample.  Therefore, we define a Gaussian distribution as the following 4-tuple:--\begin{code}--data Gaussian datapoint = Gaussian-    { n  :: {-# UNPACK #-} !Int     -- ^ The number of samples  trained on-    , m1 :: !datapoint              -- ^ The mean (first moment) of the trained distribution-    , m2 :: !datapoint              -- ^ The variance (second moment) of the trained distribution times (n-1)-    , dc :: {-# UNPACK #-} !Int     -- ^ The number of \"dummy points\" that have been added to the distribution.  Required for numerical stability reasons.-    } -    deriving (Show,Read)--{--data Gaussian = Gaussian-    { n   :: {- UNPACK -} !Int         -- the number of samples seen-    , m1  :: {- UNPACK -} !Double      -- the first moment (mean)-    , m2  :: {- UNPACK -} !Double      -- the second moment (pop. variance) times n-    , dc  :: {- UNPACK -} !Int         -- the number of "dummy" points that have been added to the Gaussian-    }-    deriving (Show,Read)--}--\end{code}--There are several important notes to make about this data structure:-\begin{enumerate}-\item Normally a Gaussian distribution would have non-negative values for both the standard deviation and the number of samples.  We make no such assumption.  The normal operations we perform on a distribution (e.g. the pdf) will only be well defined in these circumstances, but the algebraic properties of the distribution require that all parameters can be negative.--\item We have introduced a fourth variable called |dc|.  This will be necessary to track how many "dummy datapoints" we have added to the Gaussian for numerical stability reasons.  See section \ref{sec:dummy}.--%\item We do not specify the base data type to store our estimates of the parameters.  Normally, we will use the |Double| type to do this, but if we have millions of estimated Gaussians, we may want to use the |Float| type instead.  Similarly, if we require extra precision we may use the |LogFloat| data type at the expense of performance.-\end{enumerate}--For convenience, we also define the following accessor functions to the sample mean and sample variance:-\begin{code}--mean :: (Floating datatype) => Gaussian datatype -> datatype-mean (Gaussian n m1 m2 dc) = m1- -var :: (Floating datatype) => Gaussian datatype -> datatype-var (Gaussian n m1 m2 dc) = m2/(fi $ n-1)--stddev :: (Floating  datapoint, Ord datapoint) => (Gaussian datapoint) -> datapoint-stddev = sqrt . var--\end{code}--Finally, we make our Gaussian type capable of being unboxed:-\begin{code}---- derivingUnbox "Gaussian"---     [d| instance Unbox' (Gaussian) (Int, Double, Double, Int) |]---     [| \ (Gaussian n m1 m2 f) -> (n,m1,m2,f) |]---     [| \ (n,m1,m2,f) -> (Gaussian n m1 m2 f) |]---- derivingUnbox "Gaussian"---     [d| instance (U.Unbox a) => Unbox' (Gaussian a) (Int, a, a, Int) |]---     [| \ (Gaussian n m1 m2 dc) -> (n,m1,m2,dc) |]---     [| \ (n,m1,m2,dc) -> (Gaussian n m1 m2 dc) |]--derivingUnbox "Gaussian"-    [t| (U.Unbox a) => (Gaussian a) -> (Int, a, a, Int) |]-    [| \ (Gaussian n m1 m2 dc) -> (n,m1,m2,dc) |]-    [| \ (n,m1,m2,dc) -> (Gaussian n m1 m2 dc) |]---- instance NFData Gaussian where-instance (NFData datapoint) => NFData (Gaussian datapoint) where-    rnf (Gaussian n m1 m2 dc) = seq (rnf m1) $ seq (rnf m2) $ seq (rnf dc) (rnf n)--\end{code}--\section{Algebra}--We are following the homomorphic learning framework.  In this section, we will concern ourselves with algebraic manipulations of fully trained Gaussian models.  In particular, we will see how to convert them into other fully trained Gaussian models.--%We will do this by working backwards from a known batch trainer for the Gaussian distribution.  Knuth presents the %following recurrence relations in pg 232 of Vol2 AoCP:-%\begin{align*}-%m1_k &= m1_{k-1}+(x_k-m1_{k-1})/k\\-%m2_k &= m2_{k-1}+(x_k-m1_{k-1})/(x_k-m1_k)-%\end{align*}--\subsection{Semigroup}-We want to construct the semigroup operation for our Gaussian distribution so that our batch trainer will be a semigroup homomorphism.  That is, we want the property:-\begin{spec}-(train xs)<>(train ys) = train (xs++ys)-\end{spec}-To do this, we must construct appropriate definitions of |n'|, |m1'|, and |m2'| below:-\begin{spec}-(Gaussian na m1a m2a) <> (Gaussian nb m1b m2b) = Gaussian n' m1' m2'-\end{spec}-This is a somewhat long but straightforward process of applying the definitions in Section \ref{sec:defn} and then performing symbol manipulation.  We start by calculating the number of samples |n'|.  The number of samples in the resulting distribution is simply the sum of the number of samples in both the input distributions:-\begin{align}-\label{eq:n'}-n' = n_a+n_b-\end{align}--Second, we calculate the new average |m1'|.  We start with the definition that:-\begin{align}-m_1' &= \frac{1}{n'}\sum_{i=1}^{n'} x_i-\end{align}-Then we split the summation according to whether the input element $x_i$ was from the left Gaussian a or right Gaussian b, substitute equation \ref{eq:}, and rearrange:-\begin{align}-m_1' &= \frac{1}{n'}\left(\sum_{i=1}^{n_a} x_{ia} + \sum_{i=1}^{n_b} x_{ib}\right)\\-&= \frac{1}{n'}\left(n_a m_{1a} + n_b m_{1b}\right)\\-&= m_{1a}\frac{n_a}{n'} + m_{1b}\frac{n_b}{n'}-\end{align}-Notice that equation \ref{eq:m1'} is simply the weighted average of the two means.  This makes intuitive sense.  But there is a slight problem with this definition: When implemented on a computer with floating point arithmetic, we will get |Infinity| whenever |n'| is 0.  -We will solve this problem by adding a ``dummy'' element into the Gaussian whenever |n'| would be zero.  This prevents the division by zero and maintains our associativity.  See section \ref{dummy} for more info. --Finally, we must calculate the new variance |m2'|.  We start with the definition that:-\begin{align}-m_2' &= \sum_{i=1}^{n'}(x_i - m_1')^2-\end{align}-Then, we substitute with equation \ref{eq}, split the summations, and substitute with equation \ref{}:-\begin{align}-%m_2' &= \sum_{i=1}^{n'}(x_i^2 - 2x_im_1' + m_1'^2)\\-m_2' &= \sum_{i=1}^{n'}(x_i^2 - m_1'^2)\\-&= \sum_{i=1}^{n_a}(x_{ia}^2) + \sum_{i=1}^{n_b}(x_{ib}^2) - n'm_1'^2\\-\label{eq:m2'}-&= m_{2a} + n_a m_{1a}^2 + m_{2b} + n_b m_{1b}^2 - n' m_1'^2-%&= m_{2a} + \frac{m_{1a}^2}{n_a} + m_{2b} + \frac{m_{2b}^2}{n_b} - \frac{m_1'^2}{n'}-%m_2' = m_{2a}+m_{2b}+\frac{n_an_b}{n'}(m_{1a}-m_{2a})^2-\end{align}-We combine equations \ref{eq:n'}, \ref{eq:m1'} and \ref{eq:m2'} to get the semigroup instance for the |Gaussian| data type:--\subsubsection{Proof of Validity}--For this to form a valid semigroup instance, we must satisfy the associativity law by showing that for any three |g_a, g_b, g_c :: Gaussian|:-\begin{align}-(g_a \diamond g_b) \diamond g_c  = g_a \diamond (g_b \diamond g_c)-\end{align}--We do this by showing that each field in the left side of the equation ($g_l$) equals the corresponding field in the right side ($g_r$).  We can easily see this for |n| because:-$$-n_l = (n_a + n_b) + n_c = n_a + (n_b + n_c) = n_r-$$--%For |m2|, we have:-%\begin{align*}-%m_{2l} &= m_{2a} + n_a m_{1a}^2 + m_{2b} + n_b m_{2b}^2 - n' m_1'^2-%\end{align*}--%\begin{align*}-%%((m_{1a}, m_{2a}, n_a) \diamond (m_{1b}, m_{2b}, n_b)) \diamond (m_{1c}, m_{2c}, n_c)-%\end{align*}--\todo{Finish proof}---\subsubsection{The Dummy Point}-\label{dummy}-Then, when we add another Gaussian, taking the non-adjusted |n'| away from zero, we can remove the dummy element to restore the Gaussian to its original condition.  We actually depend on the group inverse to perform this operation.--We define an equivalence relation over the Gaussians that takes this into account:-\begin{align}-g_a \sim g_b \iff g_a \diamond d^i = g_b \diamond d^j-\end{align}--\subsubsection{Performance Notes}-For some reason, the commented out line results in |(<>)| an almost 10% decrease in speed.  None of the other lines cause any performance hit at all.  The commented line is not actually necessary for consistency, it's just a nicety.--The commented out line in |merge| also results in a huge performance hit.  The division must eat up a serious amount of time.-% ga == inverse gb              = identity                 -- See below.--\subsubsection{The Code}-Finally we are ready to see the actual code:--\begin{code}---- instance Semigroup Gaussian where-instance (Fractional datapoint) => Semigroup (Gaussian datapoint) where-    {-# INLINE (<>) #-}-    (<>) !ga@(Gaussian na m1a m2a fa) !gb@(Gaussian nb m1b m2b fb)-        | n'==0 && fa >0                = (inverse dummy `merge` ga) `merge` gb-        | n'==0 && fb >0                = (inverse dummy `merge` gb) `merge` ga-        | n'==0                         = (dummy `merge` ga) `merge` gb-        | n' >1 && (fa >0 || fb>0)      = (ga `merge` gb) `merge` inverse dummy-        | otherwise                     = merge ga gb-            where n' = na+nb-        -{-# INLINE dummy #-}--- dummy :: Gaussian-dummy :: (Num datatype) => Gaussian datatype-dummy = Gaussian 2 2 2 1--{-# INLINE merge #-}--- merge :: Gaussian -> Gaussian -> Gaussian-merge :: (Fractional datatype) => Gaussian datatype -> Gaussian datatype -> Gaussian datatype-merge !ga@(Gaussian na m1a m2a fa) !gb@(Gaussian nb m1b m2b fb) = Gaussian n' m1' m2' f'-    where-        n'  = na+nb-        m1' = (m1a*(fi na)+m1b*(fi nb))/(fi n')-        m2' = m2a + m2b + (m1a^2)*(fi na) + (m1b^2)*(fi nb) - (m1'^2)*(fi n')-        f' = fa+fb        --\end{code}-And our dummy data point is:-\begin{code}-\end{code}-It's possible there is a numerically optimal distribution to use for best accuracy, however, this seems to work pretty well enough in practice.---\subsection{Monoid}-Now, we must explicitly define a Monoid instance for our class as:-\begin{code}---- instance Monoid Gaussian where-instance (Fractional datapoint) => Monoid (Gaussian datapoint) where-    {-# INLINE mempty #-}-    mempty = Gaussian 0 0 0 0-    {-# INLINE mappend #-}-    mappend = (<>)--\end{code}-We will show that this is a valid left identity (the right identity is similar) by showing:-\begin{align}-%identity <> Gaussian nb m1b m2b = Gaussian n' m1' m2'-\epsilon \diamond g_b = g' \Rightarrow g_b = g'-\end{align}-This is straightforward:-\begin{align}-n' & = 0 + n_b = n_b\\-m_1' &= 0 \frac {0}{n_b} + m_{1b}\frac{n_b}{n_b} = m_{1b} \\-m_2' &= 0 + 0 + m_{2b} + n_b m_{1b}^2 - n' m_1'^2 = m_{2b} + n_b m_{1b}^2 - n_b m_{1b}^2 = m_{2b}-\end{align}--\subsection{Group}--The inverse operation for the Gaussian is defined as:-\begin{code}-  --- instance Invertible Gaussian where    -instance (Fractional datapoint) => RegularSemigroup (Gaussian datapoint) where-    {-# INLINE inverse #-}-    inverse !(Gaussian n m1 m2 dc) = Gaussian (-n) m1 (-m2) (-dc)--instance (Fractional datapoint) => Group (Gaussian datapoint)-\end{code}-We will show this is a valid left inverse (the right inverse is similar) by showing:-\begin{align}-g^{-1}\diamond g = \epsilon-\end{align}-This is straightforward:-\begin{align}-n' &= n + (-n) = 0 \\-m_1' &= 0 \text{     Only true after factoring in the dummy elements} \\-m_2' &=-\end{align}--\section{Model Training}-We use the homomorphic learning method to define the training routines for our Gaussian distribution.  That means we need to define only a single training method.  The singleton trainer is particularly simple:-\begin{code}---- | Training a Gaussian distribution takes no parameters-data GaussianParams datatype = GaussianParams--instance Model (GaussianParams datatype) (Gaussian datatype) where-    {-# INLINE getparams #-}-    getparams _ = GaussianParams-    -instance DefaultModel (GaussianParams datatype) (Gaussian datatype) where-    {-# INLINE defparams #-}-    defparams = GaussianParams--instance NFData (GaussianParams datatype) where-    rnf params = ()--instance HomTrainer (GaussianParams Double) Double (Gaussian Double) where-    train1dp' GaussianParams x = Gaussian 1 x 0 0--\end{code}--\section{Model Use}-\begin{code}-convdistr :: Gaussian Double -> N.NormalDistribution-convdistr g = N.normalDistr (mean g) (stddev g)--instance Distribution (Gaussian Double) Double Double where-    pdf g x = D.density (convdistr g) x-    cdf g x = D.cumulative (convdistr g) x-    cdfInverse g x = D.quantile (convdistr g) x-\end{code}--\section{Distribution Functions}--\subsection{Probability Density Function (PDF)}--The PDF for a Gaussian distribution is defined as:-$$-P(x) = \frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{(x-\mu)^2}{2\sigma^2}}-$$-\begin{code}-data GaussianPDF = GaussianPDF-    { prob :: Double-    , weight :: Double-    }-    -data GaussianPDFParams = GaussianPDFParams--instance Model GaussianPDFParams GaussianPDF where-    getparams _ = GaussianPDFParams-    -instance DefaultModel GaussianPDFParams GaussianPDF where-    defparams = GaussianPDFParams--instance Semigroup GaussianPDF where-    (<>) = undefined-    -instance Monoid GaussianPDF where-    mempty = GaussianPDF 0 0-    mappend = (<>)--instance HomTrainer GaussianPDFParams (Gaussian Double) GaussianPDF where-    train' GaussianPDFParams = undefined-    -\end{code}--\section{Misc}-\begin{code}--\end{code}--\section{Appendix: Symbol Manipulation}-We have defined |m2| so that:-\begin{align}-m_2 &= \sum_{i=1}^{n}(x_i - m_1)^2 \\-&= \sum_{i=1}^{n}(x_i^2 -2x_i m_1 +m_1^2)\\-&= \sum_{i=1}^{n}x_i^2 -2m_1 \sum_{i=1}^{n}x_i +\sum_{i=1}^{n}m_1^2\\-&= \sum_{i=1}^{n}x_i^2 -2m_1 n'm_1' +nm_1^2\\-\label{eq:m_2-1} &= \sum_{i=1}^{n}x_i^2 -nm_1^2\\-\label{eq:m_2-2} m_2+nm_1^2 &= \sum_{i=1}^{n}x_i^2-\end{align}--\end{document}

@@ -0,0 +1,257 @@+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE GADTs #-}++{-# LANGUAGE DataKinds #-}++module HLearn.Models.Distributions.KernelDensityEstimator+    ( KDEParams (..)+    , KDEBandwidth (..)+    , KDE (..)+    , KDE' (..)+    , genSamplePoints+    +    , Uniform (..)+    , Triangular (..)+    , Epanechnikov (..)+    , Quartic (..)+    , Triweight (..)+    , Tricube (..)+    , Gaussian (..)+    , Cosine (..)+    , KernelBox (..)+    )+    where+          +import HLearn.Algebra+import HLearn.Models.Distributions.Common+import HLearn.Models.Distributions.Categorical++import qualified Control.ConstraintKinds as CK+import Control.DeepSeq+import qualified Data.Vector.Unboxed as VU++-------------------------------------------------------------------------------+--++data KDEBandwidth prob = Constant prob | Variable (prob -> prob)++instance (Show prob) => Show (KDEBandwidth prob) where+    show (Constant x) = "Constant " ++ show x++instance (Eq prob) => Eq (KDEBandwidth prob) where+    Constant x1 == Constant x2 = x1==x2+    Variable f1 == Variable f2 = True+    _ == _ = False++instance (Ord prob) => Ord (KDEBandwidth prob) where+    Variable _ `compare` Variable _ = EQ+    Constant _ `compare` Variable _ = LT+    Variable _ `compare` Constant _ = GT+    Constant x `compare` Constant y = x `compare` y++instance (NFData prob) => NFData (KDEBandwidth prob) where+    rnf (Constant x) = rnf x+    rnf (Variable f) = seq f ()+    +{-# INLINE calcBandwidth #-}+calcBandwidth :: (KDEBandwidth prob) -> prob -> prob+calcBandwidth (Constant h) _ = h+calcBandwidth (Variable f) x = f x++data KDEParams prob = KDEParams+    { bandwidth :: KDEBandwidth prob+    , samplePoints :: VU.Vector prob -- ^ These data points must be sorted from smallest to largest+    , kernel :: KernelBox prob+    }+    deriving (Show,Eq,Ord)++instance (NFData prob) => NFData (KDEParams prob) where+    rnf kdeparams = deepseq (bandwidth kdeparams)+                  $ deepseq (kernel kdeparams)+                  $ seq (samplePoints kdeparams)+                  $ ()++genSamplePoints :: (Fractional prob, VU.Unbox prob) => Int -> Int -> Int -> VU.Vector prob+genSamplePoints min max samples = VU.fromList $ map (\i -> (fromIntegral min) + (fromIntegral i)*step) [0..samples]+    where+        step = (fromIntegral $ max-min)/(fromIntegral $ samples)++-- paramsFromData :: (Floating prob, Ord prob, Enum prob, VU.Unbox prob) => [prob] -> KDEParams prob+-- paramsFromData xs = KDEParams+--     { bandwidth = ((4*sigma^^5)/(3*n))**(1/5)+--     , samplePoints = VU.fromList samplePointsL+--     , kernel = KernelBox Epanechnikov+--     }+--     where+--         size = 100+--         margin = 10+--         step = ((maximum xs)-(minimum xs))/size+--         samplePointsL = map (\i -> (minimum xs) + i*step) [(-margin)..(size+margin)]+--         n = fromIntegral $ length xs+--         sigma = 1++data KDE' prob = KDE'+    { params :: KDEParams prob+    , n :: prob+    , sampleVals :: VU.Vector prob+    }+    deriving (Show,Eq,Ord)++instance (NFData prob) => NFData (KDE' prob) where+    rnf kde = deepseq (params kde)+            $ deepseq (n kde)+            $ seq (sampleVals kde)+            $ ()++type KDE prob = RegSG2Group (KDE' prob)++-------------------------------------------------------------------------------+-- Algebra++instance (Eq prob, Num prob, VU.Unbox prob) => Semigroup (KDE' prob) where+    kde1 <> kde2 = if (params kde1) /= (params kde2)+        then error "KDE.(<>): different params"+        else kde1+            { n = (n kde1) + (n kde2)+            , sampleVals = VU.zipWith (+) (sampleVals kde1) (sampleVals kde2)+            }++instance (Eq prob, Num prob, VU.Unbox prob) => RegularSemigroup (KDE' prob) where+    inverse kde = kde+        { n = negate $ n kde+        , sampleVals = VU.map negate $ sampleVals kde+        }++instance (Num prob, VU.Unbox prob) => LeftOperator prob (KDE' prob) where+    p .* kde = kde+        { n = p * (n kde)+        , sampleVals = VU.map (*p) (sampleVals kde)+        }++instance (Num prob, VU.Unbox prob) => RightOperator prob (KDE' prob) where+    (*.) = flip (.*)+    +-------------------------------------------------------------------------------+-- Training+    +instance (Eq prob, Num prob, VU.Unbox prob) => Model (KDEParams prob) (KDE prob) where+    getparams (SGJust kde) = params kde++instance HomTrainer (KDEParams Double) Int (KDE Double) where+    train1dp' params dp = train1dp' params (fromIntegral dp :: Double)++instance (Eq prob, Fractional prob, VU.Unbox prob) => HomTrainer (KDEParams prob) prob (KDE prob) where+    train1dp' params dp = SGJust $ KDE'+        { params = params+        , n = 1+        , sampleVals = VU.map (\x -> k ((x-dp)/(h x))) (samplePoints params)+        }+        where+            k u = (evalkernel (kernel params) u)/(h u)+            h   = calcBandwidth (bandwidth params)++-------------------------------------------------------------------------------+-- Distribution+    +instance (Ord prob, Fractional prob, VU.Unbox prob) => Distribution (KDE prob) prob prob where+    pdf (SGJust kde) dp +        | dp <= (samplePoints $ params kde) VU.! 0 = 0 -- (sampleVals kde) VU.! 0+        | dp >= (samplePoints $ params kde) VU.! l = 0 -- (sampleVals kde) VU.! l+        | otherwise = (y2-y1)/(x2-x1)*(dp-x1)+y1+        where+            index = binsearch (samplePoints $ params kde) dp+            x1 = (samplePoints $ params kde) VU.! (index-1)+            x2 = (samplePoints $ params kde) VU.! (index)+            y1 = ((sampleVals kde) VU.! (index-1)) / (n kde)+            y2 = ((sampleVals kde) VU.! (index  )) / (n kde)+            l = (VU.length $ samplePoints $ params kde)-1++binsearch :: (Ord a, VU.Unbox a) => VU.Vector a -> a -> Int+binsearch vec dp = go 0 (VU.length vec-1)+    where +        go low high+            | low==high = low+            | dp <= (vec VU.! mid) = go low mid+            | dp >  (vec VU.! mid) = go (mid+1) high+            where +                mid = floor $ (fromIntegral $ low+high)/2++-------------------------------------------------------------------------------+-- Morphisms++instance Morphism (Categorical Int Double) (KDEParams Double) (KDE Double) where+    cat $> kdeparams = train' kdeparams $ CK.fmap (fromIntegral :: Int -> Double) (cat $> FreeModParams)++-------------------------------------------------------------------------------+-- Kernels++-- | This list of kernels is take from wikipedia's: https://en.wikipedia.org/wiki/Uniform_kernel#Kernel_functions_in_common_use+class Kernel kernel num where+    evalkernel :: kernel -> num -> num++data KernelBox num where KernelBox :: (Kernel kernel num, Show kernel) => kernel -> KernelBox num++instance Kernel (KernelBox num) num where+    evalkernel (KernelBox k) p = evalkernel k p+instance Show (KernelBox num) where+    show (KernelBox k) = "KB "++show k+instance Eq (KernelBox num) where+    KernelBox k1 == KernelBox k2 = (show k1) == (show k2)+instance Ord (KernelBox num) where+    _ `compare` _ = EQ+instance NFData (KernelBox num) where+    rnf (KernelBox num)= seq num ()+    +data Uniform = Uniform deriving (Read,Show)+instance (Fractional num, Ord num) => Kernel Uniform num where+    evalkernel Uniform u = if abs u < 1+        then 1/2+        else 0++data Triangular = Triangular deriving (Read,Show)+instance (Fractional num, Ord num) => Kernel Triangular num where+    evalkernel Triangular u = if abs u<1+        then 1-abs u+        else 0+        +data Epanechnikov = Epanechnikov deriving (Read,Show)+instance (Fractional num, Ord num) => Kernel Epanechnikov num where+    evalkernel Epanechnikov u = if abs u<1+        then (3/4)*(1-u^^2)+        else 0++data Quartic = Quartic deriving (Read,Show)+instance (Fractional num, Ord num) => Kernel Quartic num where+    evalkernel Quartic u = if abs u<1+        then (15/16)*(1-u^^2)^^2+        else 0+        +data Triweight = Triweight deriving (Read,Show)+instance (Fractional num, Ord num) => Kernel Triweight num where+    evalkernel Triweight u = if abs u<1+        then (35/32)*(1-u^^2)^^3+        else 0++data Tricube = Tricube deriving (Read,Show)+instance (Fractional num, Ord num) => Kernel Tricube num where+    evalkernel Tricube u = if abs u<1+        then (70/81)*(1-u^^3)^^3+        else 0+        +data Cosine = Cosine deriving (Read,Show)+instance (Floating num, Ord num) => Kernel Cosine num where+    evalkernel Cosine u = if abs u<1+        then (pi/4)*(cos $ (pi/2)*u)+        else 0+        +data Gaussian = Gaussian deriving (Read,Show)+instance (Floating num, Ord num) => Kernel Gaussian num where+    evalkernel Gaussian u = (1/(2*pi))*(exp $ (-1/2)*u^^2)

@@ -0,0 +1,178 @@+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE FunctionalDependencies #-}++{-# LANGUAGE DataKinds #-}+{-# LANGUAGE KindSignatures #-}++{-# LANGUAGE DatatypeContexts #-}+{-# LANGUAGE NoMonomorphismRestriction #-}++module HLearn.Models.Distributions.Moments+    ( MomentsParams (..)+    , Moments (..)+    )+    where+          +import GHC.TypeLits+import qualified Data.Vector.Unboxed as VU++import HLearn.Algebra+import HLearn.Models.Distributions.Common++-------------------------------------------------------------------------------+-- Moments++data MomentsParams = MomentsParams++newtype Moments prob (n::Nat) = Moments (VU.Vector prob)+    deriving (Read,Show)++instance (Eq prob, VU.Unbox prob) => Eq (Moments prob n) where+    (Moments v1) == (Moments v2) = v1 == v2++instance (Ord prob, VU.Unbox prob) => Ord (Moments prob n) where+    compare (Moments v1) (Moments v2) = compare v1 v2++-------------------------------------------------------------------------------+-- Algebra++instance (Num prob, VU.Unbox prob) => Semigroup (Moments prob n) where+    (Moments ma) <> (Moments mb) = Moments $ VU.zipWith (+) ma mb+    +instance (Num prob, VU.Unbox prob, SingI n) => Monoid (Moments prob n) where+    mappend = (<>)+    mempty = Moments $ VU.replicate (n+1) 0+        where n=fromIntegral $ fromSing (sing :: Sing n)+    +instance (Num prob, VU.Unbox prob) => RegularSemigroup (Moments prob n) where+    inverse (Moments m) = Moments $ VU.map negate m+    +instance (Fractional prob, VU.Unbox prob) => LeftOperator prob (Moments prob n) where+    p .* (Moments vec) = Moments $ VU.map (*p) vec++instance (Fractional prob, VU.Unbox prob) => RightOperator prob (Moments prob n) where+    (*.) = flip (.*)+    +-------------------------------------------------------------------------------+-- Training+    +instance Model MomentsParams (Moments prob n) where+    getparams _ = MomentsParams++instance DefaultModel MomentsParams (Moments prob n) where+    defparams = MomentsParams++instance (VU.Unbox prob, Fractional prob, SingI n) => HomTrainer MomentsParams prob (Moments prob n) where+    train1dp' _ dp = Moments $ VU.fromList [dp^^i | i <- [0..n]]+        where n=fromIntegral $ fromSing (sing :: Sing n)+              +-------------------------------------------------------------------------------+--++data MomentsConverterParams params = MomentsConverterParams params++data MomentsConverter prob (n::Nat) dist = MomentsConverter+    { moments       :: !(Moments prob n)+    , dist          :: dist+    }+    deriving (Read,Show)++instance +    ( DefaultMorphism (Moments prob n) params dist+    , Num prob, VU.Unbox prob+    ) => Semigroup (MomentsConverter prob n dist) where+    mc1 <> mc2 = MomentsConverter m (morph m)+        where+            m = (moments mc1) <> (moments mc2)++-------------------------------------------------------------------------------+--++data BetaParams = BetaParams++data Beta prob = Beta+    { alpha :: prob+    , beta :: prob+    }+    deriving (Read,Show)++instance (VU.Unbox prob, Fractional prob) => Morphism (Moments prob 2) (BetaParams) (Beta prob) where+    (Moments v) $> BetaParams = Beta+        { alpha = alpha+        , beta  = beta+        }+        where+            alpha = (k-1)/(l*(1+k))+            beta  = k*alpha+            +            k = (mean-1)/mean+            l = ((k+1)^^2)*var+            +            mean = 1+            var = 1++-------------------------------------------------------------------------------+--+++data NormalParams = NormalParams+              +data Normal prob = Normal+    { n :: prob+    , mean :: prob+    , stddev :: prob+    }+    deriving (Read,Show)+    +{-instance Semigroup (Normal prob)+instance Monoid (Normal prob)+    +instance Model NormalParams (Normal prob) where+    getparams normal = NormalParams+    +instance DefaultModel NormalParams (Normal prob) where+    defparams = NormalParams-}+    +instance (VU.Unbox prob, Fractional prob) => Morphism (Moments prob 2) NormalParams (Normal prob) where+-- instance Morphism (Moments prob 2) NormalParams (Normal prob) where+    (Moments v) $> NormalParams = Normal+        { n         = m0+        , mean      = m1 / m0+        , stddev    = (1/(m0-1))*m2-(m0/(m0-1))*(m1/m0)^^2+        }+        where+            m0 = v VU.! 0+            m1 = v VU.! 1+            m2 = v VU.! 2+            +instance (VU.Unbox prob, Fractional prob) => Morphism (Normal prob) MomentsParams (Moments prob 2) where+    (Normal n mean stddev) $> MomentsParams = Moments $ VU.fromList+        [ n+        , mean * n+        , (stddev+(n/(n-1))*(mean)^^2)*(n-1)+        ]++foo = ((train' MomentsParams [1,2,3::Double] :: Moments Double 2)+    $> NormalParams :: Normal Double)++{-foo' = train' +    ( (MomentsParams :: MomentsParams)+    ) [1,2,3]+    $> NormalParams :: Normal Double-}+    +-- foo2 = (train' +--     ( ((NormalParams :: NormalParams)+--     :. (MomentsParams :: MomentsParams))+--         :: (MorphismComposition+--                           [Double]+--                           (MomentsParams)+--                           (Moments Double 2)+--                           (NormalParams)+--                           (Normal Double))+--     ) [1,2,3::Double]) :: Normal Double

@@ -13,12 +13,13 @@ import qualified Control.ConstraintKinds as CK  size = 10^8-main = defaultMain +main = putStrLn "done."+{-main = defaultMain      [ bench "HLearn-Gaussian" $ nf ((train :: VU.Vector Double -> Gaussian Double)) (VU.enumFromN (0::Double) size)     , bench "HLearn-Gaussian-Parallel" $ whnf (parallel $ (train :: VU.Vector Double -> Gaussian Double)) (VU.enumFromN (0::Double) size) --     , bench "HLearn-Gaussian-List" $ nf (train :: [Double] -> Gaussian Double) [0..fromIntegral size]     , bench "statistics-Gaussian" $ whnf (normalFromSample . VU.enumFromN 0) (size)-    ]+    ]-}          {-bench "batch train [] 1e6" $ nf ((batch train) GaussianParams) [0..1e6::Double] --     , bench "batch train V 1e6"  $ nf ((batch train) GaussianParams) (V.enumFromN (0::Double) (10^6))

@@ -18,6 +18,7 @@ size=10^8 main = do --     print $ (train (concat $ replicate (10^6) ["a","b"]) :: Categorical String Double)-    print $ (train [1..fromIntegral size::Double] :: Gaussian Double)-    print $ (train (VU.enumFromN (0::Double) size) :: Gaussian Double)+--     print $ (train [1..fromIntegral size::Double] :: Gaussian Double)+--     print $ (train (VU.enumFromN (0::Double) size) :: Gaussian Double) --     print $ normalFromSample $ VU.enumFromN 0 size+    putStrLn "done."