diff --git a/HLearn-distributions.cabal b/HLearn-distributions.cabal --- a/HLearn-distributions.cabal +++ b/HLearn-distributions.cabal @@ -1,9 +1,9 @@ Name: HLearn-distributions -Version: 0.2.2.1 +Version: 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. +Description: This module is used to estimate statistical distributions from data. It is based on the algebraic properties of the "HomTrainer" type class from the HLearn-algebra package. Category: Data Mining, Machine Learning, Statistics -License: GPL +License: BSD3 --License-file: LICENSE Author: Mike izbicki Maintainer: mike@izbicki.me @@ -12,54 +12,13 @@ homepage: http://github.com/mikeizbicki/HLearn/ bug-reports: http://github.com/mikeizbicki/HLearn/issues -Executable HLearn-Distributions-Criterion - Main-is: src/examples/Criterion.hs - Build-Depends: - HLearn-algebra >= 0.0.1, - ConstraintKinds >= 0.0.1, - HLearn-distributions >= 0.0.1, - - base >= 3 && < 5, - criterion >= 0.6.1.1, - vector, --- logfloat , - statistics - ghc-options: - -threaded - -rtsopts - -O2 - -funbox-strict-fields - -- -prof - -- -fllvm - -Executable HLearn-Distributions-SpaceTests - Main-is: src/examples/SpaceTests.hs - Build-Depends: - HLearn-algebra >= 0.0.1, - --HLearn-algebra , - ConstraintKinds, - - HLearn-distributions , - base >= 3 && < 5, - criterion >= 0.6.1.1, - vector, - logfloat , - statistics - ghc-options: - -threaded - -rtsopts - -O2 - -funbox-strict-fields - --enable-executable-profiling - -- -prof - -- -fllvm - Library Build-Depends: HLearn-algebra >= 0.1.2, ConstraintKinds >= 0.0.1, base >= 3 && < 5, + template-haskell, deepseq >= 1.3.0.1, list-extras >= 0.4.1, containers >= 0.5, @@ -68,8 +27,11 @@ QuickCheck >= 2.5.1, vector >= 0.9, vector-th-unbox >= 0.2, + graphviz >= 2999.16, + hmatrix >= 0.14, -- are these really necessary? + array >= 0.4.0, process >= 1.1.0.2, MonadRandom >= 0.1.6, math-functions >= 0.1.1, @@ -85,17 +47,26 @@ -O2 -- -fllvm Exposed-modules: - HLearn.Gnuplot.Distributions HLearn.Models.Distributions HLearn.Models.Distributions.Common - HLearn.Models.Distributions.Categorical - HLearn.Models.Distributions.KernelDensityEstimator - HLearn.Models.Distributions.KernelDensityEstimator.Kernels - HLearn.Models.Distributions.Moments - HLearn.Models.Distributions.Multivariate - --HLearn.Models.Distributions.Normal - HLearn.Models.Distributions.Gaussian - --HLearn.Models.Distributions.GaussianOld - --HLearn.Models.Distributions.GaussianOld2 - --HLearn.Models.Distributions.Poisson + HLearn.Models.Distributions.Visualization.Gnuplot + HLearn.Models.Distributions.Visualization.Graphviz + HLearn.Models.Distributions.Univariate.Binomial + HLearn.Models.Distributions.Univariate.Categorical + HLearn.Models.Distributions.Univariate.Exponential + HLearn.Models.Distributions.Univariate.Geometric + HLearn.Models.Distributions.Univariate.LogNormal + HLearn.Models.Distributions.Univariate.Normal + HLearn.Models.Distributions.Univariate.Poisson + --HLearn.Models.Distributions.Univariate.Uniform + HLearn.Models.Distributions.Univariate.Internal.MissingData + HLearn.Models.Distributions.Univariate.Internal.Moments + HLearn.Models.Distributions.Multivariate.Interface + HLearn.Models.Distributions.Multivariate.MultiNormal + HLearn.Models.Distributions.Multivariate.Internal.CatContainer + HLearn.Models.Distributions.Multivariate.Internal.Container + HLearn.Models.Distributions.Multivariate.Internal.Ignore + HLearn.Models.Distributions.Multivariate.Internal.Marginalization + HLearn.Models.Distributions.Multivariate.Internal.TypeLens + HLearn.Models.Distributions.Multivariate.Internal.Unital diff --git a/src/HLearn/Gnuplot/Distributions.hs b/src/HLearn/Gnuplot/Distributions.hs deleted file mode 100644 --- a/src/HLearn/Gnuplot/Distributions.hs +++ /dev/null @@ -1,133 +0,0 @@ -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE TypeSynonymInstances #-} - --- | This file contains the functions for plotting distributions using Gnuplot. The interface is still under heavy construction, so it's not very well documented as of yet. - -module HLearn.Gnuplot.Distributions - where - -import HLearn.Algebra -import HLearn.Models.Distributions - -import qualified Data.Map as Map -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 - } - -plotFile :: String -> PlotParams -plotFile str = PlotParams - { dataFile = str++".dat" - , gnuFile = str++".gnu" - , picFile = str++".ps" - } - --- | provided due to backwards compatibility with the nuclear weapons blog post. -genPlotParams :: String -> a -> PlotParams -genPlotParams str a = plotFile str - -defPlotParams = PlotParams - { dataFile = "hlearn-distributions.dat" - , gnuFile = "hlearn-distributions.gnu" - , picFile = "hlearn-distributions.ps" - } - -class PlottableDistribution dist where - plotDistribution :: PlotParams -> dist -> IO () - plotDistribution params dist = do - -- Create data file - putStrLn "Creating data file..." - datah <- openFile (dataFile params) WriteMode - hPutStrLn datah $ plotdata dist - hClose datah - - -- Create gnuplot file - putStrLn "Plotting data" - gnuh <- openFile (gnuFile params) WriteMode - hPutStrLn gnuh $ gnuplot params dist - hClose gnuh - - -- Run gnuplot, creating picture - system $ "gnuplot "++(gnuFile params) - putStrLn "done." - return () - - plotdata :: dist -> String - gnuplot :: PlotParams -> dist -> String - -instance PlottableDistribution (KDE Double) where - plotdata dist = mconcat [show (x::Double) ++ " " ++ show (pdf dist x::Double) | x <- plotPoints] - where - plotPoints = VU.toList $ samplePoints $ getparams dist - - gnuplot params dist - = "set terminal postscript \"Times-Roman\" 25" - ++ "set output \"" ++ (picFile params) ++ "\"" - ++ "unset xtics; unset ytics; unset key" - ++ "set border 0; set xzeroaxis lt 1; set yzeroaxis lt 1" - ++ "plot '"++(dataFile params)++"' using 1:2 lt 1 lw 4 lc rgb '#ccccff' with filledcurves, \\" - ++ " '"++(dataFile params)++"' using 1:2 lt 1 lw 4 lc rgb '#0000ff' with lines" - -instance (Show label, Show prob, Num prob, Ord prob) => PlottableDistribution (Categorical label prob) where - plotdata dist@(Categorical pdf) = concat $ do - (label,prob) <- Map.toList pdf - return $ show label++" "++show prob++"\n" - - gnuplot params (Categorical pdf) - = "set terminal postscript \"Times-Roman\" 25; set output \""++ (picFile params) ++ "\"\n" - ++"set tics scale 0; set xtics nomirror; unset ytics; unset key\n" - ++"set border 2; set xzeroaxis lt 1\n" - ++"set ylabel \"Probability\"\n" - ++yrange - ++"set style fill border -1\n" - ++"set style data histogram; set style histogram cluster gap 1\n" - ++"plot '"++(dataFile params)++"' using 2:xticlabels(1)\n" - where - positiveSamples = or $ map (\(k,v) -> v>0) $ Map.toList pdf - negativeSamples = or $ map (\(k,v) -> v<0) $ Map.toList pdf - yrange = if positiveSamples && negativeSamples - then "" - else if positiveSamples - then "set yrange [0:]\n" - else "set yrange [:0]\n" -{- = "set terminal postscript \"Times-Roman\" 25 \n" - ++ "set output \"" ++ (picFile params) ++ "\" \n" - ++ "unset xtics; unset ytics; unset key\n" - ++ "set border 0; set xzeroaxis lt 1; set yzeroaxis lt 1\n" - ++ "set style data histogram \n" - ++ "plot '"++(dataFile params)++"' using 1:2 \n"-} - --- 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 () diff --git a/src/HLearn/Models/Distributions.hs b/src/HLearn/Models/Distributions.hs --- a/src/HLearn/Models/Distributions.hs +++ b/src/HLearn/Models/Distributions.hs @@ -7,17 +7,36 @@ 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.Moments - , module HLearn.Models.Distributions.Multivariate + , module HLearn.Models.Distributions.Visualization.Gnuplot + , module HLearn.Models.Distributions.Visualization.Graphviz + , module HLearn.Models.Distributions.Univariate.Binomial + , module HLearn.Models.Distributions.Univariate.Categorical + , module HLearn.Models.Distributions.Univariate.Exponential + , module HLearn.Models.Distributions.Univariate.Geometric + , module HLearn.Models.Distributions.Univariate.LogNormal + , module HLearn.Models.Distributions.Univariate.Normal +-- , module HLearn.Models.Distributions.Univariate.Uniform + , module HLearn.Models.Distributions.Univariate.Poisson + , module HLearn.Models.Distributions.Univariate.Internal.MissingData +-- , module HLearn.Models.Distributions.KernelDensityEstimator + , module HLearn.Models.Distributions.Multivariate.Interface + , module HLearn.Models.Distributions.Multivariate.MultiNormal + , module HLearn.Models.Distributions.Multivariate.Internal.TypeLens ) 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.Moments -import HLearn.Models.Distributions.Multivariate +import HLearn.Models.Distributions.Visualization.Gnuplot +import HLearn.Models.Distributions.Visualization.Graphviz +import HLearn.Models.Distributions.Univariate.Binomial +import HLearn.Models.Distributions.Univariate.Categorical +import HLearn.Models.Distributions.Univariate.Exponential +import HLearn.Models.Distributions.Univariate.Geometric +import HLearn.Models.Distributions.Univariate.LogNormal +import HLearn.Models.Distributions.Univariate.Normal +-- import HLearn.Models.Distributions.Univariate.Uniform +import HLearn.Models.Distributions.Univariate.Poisson +import HLearn.Models.Distributions.Univariate.Internal.MissingData +import HLearn.Models.Distributions.Multivariate.Interface +import HLearn.Models.Distributions.Multivariate.MultiNormal +import HLearn.Models.Distributions.Multivariate.Internal.TypeLens diff --git a/src/HLearn/Models/Distributions/Categorical.hs b/src/HLearn/Models/Distributions/Categorical.hs deleted file mode 100644 --- a/src/HLearn/Models/Distributions/Categorical.hs +++ /dev/null @@ -1,147 +0,0 @@ -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE UndecidableInstances #-} -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE BangPatterns #-} - - --- | The categorical distribution is used for discrete data. It is also sometimes called the discrete distribution or the multinomial distribution. For more, see the wikipedia entry: -module HLearn.Models.Distributions.Categorical - ( - -- * Data types - Categorical (..) - , CategoricalParams(..) - - -- * Helper functions - , dist2list - , mostLikely - ) - where - -import Control.DeepSeq -import Control.Monad.Random -import Data.List -import Data.List.Extras -import Debug.Trace - -import qualified Data.Map.Strict as Map -import qualified Data.Foldable as F - -import HLearn.Algebra -import HLearn.Models.Distributions.Common - -------------------------------------------------------------------------------- --- CategoricalParams - --- | The Categorical distribution takes no parameters -data CategoricalParams = CategoricalParams - deriving (Read,Show,Eq,Ord) - -instance NFData CategoricalParams where - rnf x = () - -instance Model CategoricalParams (Categorical label probtype) where - getparams model = CategoricalParams - -instance DefaultModel CategoricalParams (Categorical label probtype) where --- instance DefaultModel CategoricalParams (Categorical Int Double) where - defparams = CategoricalParams - -------------------------------------------------------------------------------- --- Categorical - -data Categorical sampletype probtype = Categorical - { pdfmap :: !(Map.Map sampletype probtype) - } - deriving (Show,Read,Eq,Ord) - - -instance (NFData sampletype, NFData probtype) => NFData (Categorical sampletype probtype) where - rnf d = rnf $ pdfmap d - -------------------------------------------------------------------------------- --- Training - -instance (Ord label, Num probtype) => HomTrainer CategoricalParams label (Categorical label probtype) where - train1dp' params dp = Categorical $ Map.singleton dp 1 - -------------------------------------------------------------------------------- --- Distribution - -instance (Ord label, Ord prob, Fractional prob, Random prob) => Distribution (Categorical label prob) label prob where - - {-# INLINE pdf #-} - pdf dist label = {-0.0001+-}(val/tot) - where - val = case Map.lookup label (pdfmap dist) of - Nothing -> 0 - Just x -> x - tot = F.foldl' (+) 0 $ pdfmap dist - -{- {-# INLINE cdf #-} - cdf dist label = (Map.foldl' (+) 0 $ Map.filterWithKey (\k a -> k<=label) $ pdfmap dist) - / (Map.foldl' (+) 0 $ pdfmap dist) - - {-# INLINE cdfInverse #-} - cdfInverse dist prob = go prob pdfL - where - pdfL = map (\k -> (k,pdf dist k)) $ Map.keys $ pdfmap dist - go prob [] = fst $ last pdfL - go prob (x:xs) = if prob < snd x && prob > (snd $ head xs) - then fst x - else go (prob-snd x) xs --- cdfInverse dist prob = argmax (cdf dist) $ Map.keys $ pdfmap dist - - {-# INLINE mean #-} - mean dist = fst $ argmax snd $ Map.toList $ pdfmap dist - - {-# INLINE drawSample #-} - drawSample dist = do - x <- getRandomR (0,1) - return $ cdfInverse dist (x::prob) --} - --- | Extracts the element in the distribution with the highest probability -mostLikely :: Ord prob => Categorical label prob -> label -mostLikely dist = fst $ argmax snd $ Map.toList $ pdfmap dist - --- | Converts a distribution into a list of (sample,probability) pai -dist2list :: Categorical sampletype probtype -> [(sampletype,probtype)] -dist2list (Categorical pdfmap) = Map.toList pdfmap - -------------------------------------------------------------------------------- --- 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 - res = Map.unionWith (+) (pdfmap d1) (pdfmap d2) - -instance (Ord label, Num probtype) => RegularSemigroup (Categorical label probtype) where - inverse d1 = d1 {pdfmap=Map.map (0-) (pdfmap d1)} - -instance (Ord label, Num probtype) => Monoid (Categorical label probtype) where - mempty = Categorical Map.empty - mappend = (<>) - --- 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 diff --git a/src/HLearn/Models/Distributions/Common.hs b/src/HLearn/Models/Distributions/Common.hs --- a/src/HLearn/Models/Distributions/Common.hs +++ b/src/HLearn/Models/Distributions/Common.hs @@ -1,21 +1,106 @@ -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE FunctionalDependencies #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE EmptyDataDecls #-} +-- | This module contains the type classes for manipulating distributions. +-- +-- We use the same classes for both discrete and continuous distributions. Unfortunately, we cannot use the type classes from the 'statistics' package because we require more flexibility than they offer. + module HLearn.Models.Distributions.Common - ( Distribution(..) + ( + -- * Type classes + Probabilistic(..) + , CDF(..) + , PDF(..) + , Mean(..) + , Variance(..) + + -- * Utility functions + , nonoverlap ) where +import Data.List import HLearn.Algebra ------------------------------------------------------------------------------- -- Distribution + +-- | +class Probabilistic model where + type Probability model + +-- | Technically, every distribution has a Cumulative Distribution Function (CDF), and so this type class should be merged with the "Distribution" type class. However, I haven't had a chance to implement the CDF for most distributions yet, so this type class has been separated out. +class (Probabilistic dist) => CDF dist where +-- class CDF dist dp prob | dist -> dp, dist -> prob where + cdf :: dist -> Datapoint dist -> Probability dist + cdfInverse :: dist -> Probability dist -> Datapoint dist + +-- | Not every distribution has a Probability Density Function (PDF), however most distributions in the HLearn library do. For many applications, the PDF is much more intuitive and easier to work with than the CDF. For discrete distributions, this is often called a Probability Mass Function (PMF); however, for simplicity we use the same type class for both continuous and discrete data. +class (Probabilistic dist) => PDF dist where + pdf :: dist -> Datapoint dist -> Probability dist + +class (HomTrainer dist) => Mean dist where + mean :: dist -> Datapoint dist --- | 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 dp prob | dist -> dp, dist -> prob where - pdf :: dist -> dp -> prob --- cdf :: dist -> dp -> prob --- cdfInverse :: dist -> prob -> dp +class (Probabilistic dist) => Variance dist where + variance :: dist -> Probability dist + +-- class PDF dist dp prob | dist -> dp, dist -> prob where +-- pdf :: dist -> dp -> prob + + -- mean :: dist -> sample -- drawSample :: (RandomGen g) => dist -> Rand g sample + +-- class (Distribution dist dp prob) => DistributionOverlap dist dp prob where +-- overlap :: [dist] -> prob + +-- instance DistributionOverlap dist Double prob where +-- overlap xs = fmap (sort . (flip pdf) [-10..10]) xs + +------------------------------------------------------------------------------- +-- Continuity + +data Discrete +data Continuous + +type family Continuity t :: * +type instance Continuity Int = Discrete +type instance Continuity Integer = Discrete +type instance Continuity Char = Discrete +type instance Continuity String = Discrete + +type instance Continuity Float = Continuous +type instance Continuity Double = Continuous +type instance Continuity Rational = Continuous + +------------------------------------------------------------------------------- +-- Utilities + +-- | If you were to plot a list of distributions, nonoverlap returns the amount of area that only a single distribution covers. That is, it will be equal to number of distributions - the overlap. +-- +-- This function is used by the HomTree classifier. +nonoverlap :: + ( Enum (Probability dist), Fractional (Probability dist), Ord (Probability dist) + , PDF dist + , CDF dist + ) => [dist] -> Probability dist +nonoverlap xs = (sum scoreL)/(fromIntegral $ length scoreL) + where + scoreL = fmap (diffscore . reverse . sort . normalizeL) $ transpose $ fmap ((flip fmap) sampleL . pdf) xs + samplelen = length sampleL + sampleL = concat $ fmap sampleDist xs + sampleDist dist = fmap (cdfInverse dist . (/100)) [1..99] + +diffscore :: (Num prob) => [prob] -> prob +diffscore (x1:x2:xs) = x1-x2 + +weightedscore :: (Num prob) => [prob] -> prob +weightedscore = sum . fmap (\(i,v) -> (fromIntegral i)*v) . zip [0..] + +normalizeL :: (Fractional a) => [a] -> [a] +normalizeL xs = fmap (/tot) xs + where + tot = sum xs diff --git a/src/HLearn/Models/Distributions/Gaussian.lhs b/src/HLearn/Models/Distributions/Gaussian.lhs deleted file mode 100644 --- a/src/HLearn/Models/Distributions/Gaussian.lhs +++ /dev/null @@ -1,466 +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: . - -module HLearn.Models.Distributions.Gaussian - ( Gaussian - , GaussianParams (..) - ) - where - -import Control.Monad -import Control.Monad.Random -import Control.DeepSeq -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,Eq,Ord) - -{- -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 - deriving (Read,Show,Eq,Ord) - -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} diff --git a/src/HLearn/Models/Distributions/KernelDensityEstimator.hs b/src/HLearn/Models/Distributions/KernelDensityEstimator.hs deleted file mode 100644 --- a/src/HLearn/Models/Distributions/KernelDensityEstimator.hs +++ /dev/null @@ -1,207 +0,0 @@ -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE UndecidableInstances #-} -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE BangPatterns #-} -{-# LANGUAGE TypeOperators #-} -{-# LANGUAGE FunctionalDependencies #-} -{-# LANGUAGE RankNTypes #-} -{-# LANGUAGE GADTs #-} - -{-# LANGUAGE DataKinds #-} - --- | Kernel Density Estimation (KDE) is a generic and powerful method for estimating a probability distribution. See wikipedia for more information: -module HLearn.Models.Distributions.KernelDensityEstimator - ( - -- * Parameters - KDEParams (..) - , KDEBandwidth (..) - - -- * Data types - , KDE --- , KDE' - - -- * Easy creation of parameters - , genSamplePoints - --- , module HLearn.Models.Distributions.KernelDensityEstimator.Kernels - ) - where - -import HLearn.Algebra -import HLearn.Models.Distributions.Common -import HLearn.Models.Distributions.Categorical -import HLearn.Models.Distributions.KernelDensityEstimator.Kernels - -import qualified Control.ConstraintKinds as CK -import Control.DeepSeq -import qualified Data.Vector.Unboxed as VU - -------------------------------------------------------------------------------- --- Parameters - --- | The bandwidth is a \"free parameter\" that has a large influence on the resulting PDF. The simplest way to set the bandwidth is using a constant. However, we can also have a bandwidth that varies over the domain of the distribution. This is done using the Variable constructor. For more, see Wikipedia's entry on Variable KDE: . - -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) - $ () - - --- | Generates a vector of the positions to sample our distribution at to generate the PDF. It is intended for use with KDEParams only. -genSamplePoints :: - ( Fractional prob - , VU.Unbox prob - ) => Int -- minimum point - -> Int -- maximum point - -> Int -- number of samples - -> 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 - -------------------------------------------------------------------------------- --- KDE - --- | This is an intermediate data structure used by the KDE data type. You should NEVER use it directly because it doesn't have an empty element. It is exported because some other modules in the HLearn library need direct access to it. -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) - $ () - --- | The data structure that stores our Kernel Density Estimate. Because KDE' doesn't have an empty element, we use the RegSG2Group free structure to provide one for us. KDE inherits all the instances of KDE', plus Monoid and Group instances. -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) - diff --git a/src/HLearn/Models/Distributions/KernelDensityEstimator/Kernels.hs b/src/HLearn/Models/Distributions/KernelDensityEstimator/Kernels.hs deleted file mode 100644 --- a/src/HLearn/Models/Distributions/KernelDensityEstimator/Kernels.hs +++ /dev/null @@ -1,89 +0,0 @@ -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE GADTs #-} -{-# LANGUAGE FlexibleInstances #-} - --- | This list of kernels is take from wikipedia's: -module HLearn.Models.Distributions.KernelDensityEstimator.Kernels - ( - -- * Data types - Kernel (..) - , KernelBox (..) - - -- * Kernels - , Uniform (..) - , Triangular (..) - , Epanechnikov (..) - , Quartic (..) - , Triweight (..) - , Tricube (..) - , Gaussian (..) - , Cosine (..) - - ) - where - -import Control.DeepSeq - --- | A kernel is function in one parameter that takes a value on the x axis and spits out a "probability." We create a data object for each kernel, and a corresponding class to make things play nice with the type system. -class Kernel kernel num where - evalkernel :: kernel -> num -> num - --- | A KernelBox is a universal object for storing kernels. Whatever kernel it stores, it becomes a kernel with the same properties. -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) diff --git a/src/HLearn/Models/Distributions/Moments.hs b/src/HLearn/Models/Distributions/Moments.hs deleted file mode 100644 --- a/src/HLearn/Models/Distributions/Moments.hs +++ /dev/null @@ -1,183 +0,0 @@ -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE UndecidableInstances #-} -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE BangPatterns #-} -{-# LANGUAGE TypeOperators #-} -{-# LANGUAGE FunctionalDependencies #-} - -{-# LANGUAGE DataKinds #-} -{-# LANGUAGE KindSignatures #-} - -{-# LANGUAGE DatatypeContexts #-} -{-# LANGUAGE NoMonomorphismRestriction #-} - --- | The method of moments can be used to estimate a number of commonly used distributions. This module is still under construction as I work out the best way to handle morphisms from the Moments type to types of other distributions. For more information, see the wikipedia entry: - -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) - --- instance Distribution (Moments Double 2) Double Double where --- pdf dist dp = - -{-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 diff --git a/src/HLearn/Models/Distributions/Multivariate.hs b/src/HLearn/Models/Distributions/Multivariate.hs deleted file mode 100644 --- a/src/HLearn/Models/Distributions/Multivariate.hs +++ /dev/null @@ -1,187 +0,0 @@ -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE UndecidableInstances #-} -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE BangPatterns #-} -{-# LANGUAGE TypeOperators #-} -{-# LANGUAGE FunctionalDependencies #-} -{-# LANGUAGE ConstraintKinds #-} - -{-# LANGUAGE DataKinds #-} -{-# LANGUAGE KindSignatures #-} -{-# LANGUAGE ExistentialQuantification #-} -{-# LANGUAGE RankNTypes #-} --- {-# LANGUAGE DatatypeContexts #-} - --- | - -module HLearn.Models.Distributions.Multivariate -{- ( MultivariateParams (..) - , Multivariate - , (:::) (..) - )-} - where - -import HLearn.Algebra -import HLearn.Models.Distributions.Common -import HLearn.Models.Distributions.Categorical -import GHC.TypeLits - --- data Discrete --- data Continuous --- data Copula - --- data ParamsBox = forall (Model params dist) => params . ParamsBox params --- data DistBox - -------------------------------------------------------------------------------- --- Crazy Tuple - --- data T a b c = (a c) ::: (b c) --- deriving (Read,Show,Eq,Ord) - -data a ::: b = a ::: b - deriving (Read,Show,Eq,Ord) - -instance (Semigroup a, Semigroup b) => Semigroup (a ::: b) where - (a1 ::: b1) <> (a2 ::: b2) = (a1<>a2):::(b1<>b2) - -instance (Monoid a, Monoid b) => Monoid (a ::: b) where - mempty = mempty ::: mempty - (a1 ::: b1) `mappend` (a2 ::: b2) = (a1 `mappend` a2):::(b1 `mappend` b2) - -instance (RegularSemigroup a, RegularSemigroup b) => RegularSemigroup (a ::: b) where - inverse (a ::: b) = (inverse a):::(inverse b) - -------------------------------------------------------------------------------- --- Multivariate - -data MultivariateParams distparams {-copulaparams-} = MultivariateParams - { distparams :: distparams --- , copulaparams :: copulaparams - } - deriving (Read,Show,Eq,Ord) - -data Multivariate' distparams {-copulaparams-} distL {-copula-} = Multivariate' - { params :: MultivariateParams distparams {-copulaparams-} - , distL :: distL --- , copula :: copula - } - deriving (Read,Show,Eq,Ord) - --- data Multivariate' distparams {-copulaparams-} distL {-copula-} = Multivariate' --- { params :: MultivariateParams distparams {-copulaparams-} --- , distL :: distL --- -- , copula :: copula --- } --- deriving (Read,Show,Eq,Ord) - -data ParamsMask dist = forall params dist . (Model params dist) => ParamsMask dist - --- type Multivariate distL = forall distparams distL . (Model distparams distL) => RegSG2Group (Multivariate' distparams distL) -type Multivariate distparams distL = RegSG2Group (Multivariate' distparams distL) --- Multivariate distL = forall (Model distparams distL) => RegSG2Group (Multivariate' distparams distL) - --- data Booger distL = forall distparams distL . (Model distparams distL) => Booger (Multivariate distparams distL) - -------------------------------------------------------------------------------- --- Algebra - -instance (Semigroup distL, Eq distparams) => Semigroup (Multivariate' distparams distL) where - mv1 <> mv2 = if (params mv1)/=(params mv2) - then error "Multivariate.(<>): Adding distributions with different params" - else mv1 { distL = (distL mv1) <> (distL mv2) } - -instance (RegularSemigroup distL, Eq distparams) => RegularSemigroup (Multivariate' distparams distL) where - inverse mv = mv { distL = inverse $ distL mv } - -------------------------------------------------------------------------------- --- Training - -instance - ( Model params1 model1 - , Model params2 model2 - ) => Model (params1 ::: params2) (model1 ::: model2) where - getparams (model1 ::: model2) = (getparams model1 ::: getparams model2) - -instance - ( DefaultModel params1 model1 - , DefaultModel params2 model2 - , Model (params1 ::: params2) (model1 ::: model2) - ) => DefaultModel (params1 ::: params2) (model1 ::: model2) where - defparams = (defparams ::: defparams) - -instance (RegularSemigroup distL, Eq distparams) => Model (MultivariateParams distparams) (Multivariate distparams distL) where - getparams (SGJust mv) = params mv - -instance - ( DefaultModel params dist - , DefaultModel paramsT distT - , RegularSemigroup dist - , RegularSemigroup distT - , Eq params - , Eq paramsT - ) => DefaultModel (MultivariateParams (params ::: paramsT)) (Multivariate (params ::: paramsT) (dist ::: distT)) - where - defparams = MultivariateParams (defparams ::: defparams) - -instance - ( HomTrainer params1 dp1 model1 - , HomTrainer params2 dp2 model2 - ) => HomTrainer (params1 ::: params2) (dp1 ::: dp2) (model1 ::: model2) - where - train1dp' (params1 ::: params2) (dp1 ::: dp2) = (train1dp' params1 dp1 ::: train1dp' params2 dp2) - -instance - ( HomTrainer params dp dist - , HomTrainer paramsT dpT distT - , RegularSemigroup dist - , RegularSemigroup distT - , Eq paramsT - , Eq params - ) => HomTrainer (MultivariateParams (params ::: paramsT)) (dp ::: dpT) (Multivariate (params ::: paramsT) (dist ::: distT)) - where - train1dp' (MultivariateParams (params ::: paramsT)) (dp ::: dpT) = - SGJust $ Multivariate' (MultivariateParams (params ::: paramsT)) (train1dp' params dp ::: train1dp' paramsT dpT) - -------------------------------------------------------------------------------- --- Distribution - -instance - ( Eq params - , RegularSemigroup dist - , Distribution dist dp prob - ) => Distribution (Multivariate params dist) dp prob where - pdf (SGJust dist) = pdf (distL dist) - -instance - ( Num prob - , Distribution dist dp prob - , Distribution distT dpT prob - ) => Distribution (dist:::distT) (dp:::dpT) prob where - pdf (dist:::distT) (dp:::dpT) = (pdf dist dp)*(pdf distT dpT) - --- instance --- ( Eq params --- , Eq paramsT --- , RegularSemigroup dist --- , RegularSemigroup distT --- , Distribution dist dp prob --- ) => Distribution (Multivariate (params ::: paramsT) (dist:::distT)) (dp:::dpT) prob where --- pdf (SGJust (Multivariate' params (dist:::distT))) (dp ::: dpT) = (pdf dist dp)*(pdf () dpT) - -------------------------------------------------------------------------------- --- Testing - -test = train1dp ("food":::6:::4) - :: Multivariate (CategoricalParams ::: CategoricalParams ::: CategoricalParams) - (Categorical String Double ::: Categorical Int Double ::: Categorical Int Double) - -test2 = train [("food" ::: "dog") - ,("food" ::: "cat") - ,("drink" ::: "water") - ,("drink" ::: "piss") - ] - :: Multivariate (CategoricalParams ::: CategoricalParams) (Categorical String Double ::: Categorical String Double) - diff --git a/src/HLearn/Models/Distributions/Multivariate/Interface.hs b/src/HLearn/Models/Distributions/Multivariate/Interface.hs new file mode 100644 --- /dev/null +++ b/src/HLearn/Models/Distributions/Multivariate/Interface.hs @@ -0,0 +1,145 @@ +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE FunctionalDependencies #-} +{-# LANGUAGE PolyKinds #-} +{-# LANGUAGE StandaloneDeriving #-} + +-- | Used for Multivariate distributions + +module HLearn.Models.Distributions.Multivariate.Interface + ( + Multivariate + + -- * Type functions +-- , Ignore + , MultiCategorical (..) + , Independent (..) + , Dependent (..) + + -- * Modules + , Index (..) + , module HLearn.Models.Distributions.Multivariate.Internal.Ignore + , module HLearn.Models.Distributions.Multivariate.Internal.Marginalization + ) + where + +import Control.DeepSeq +import GHC.TypeLits + +import HLearn.Algebra +import HLearn.Models.Distributions.Common +import HLearn.Models.Distributions.Multivariate.Internal.CatContainer hiding (ds,baseparams) +import HLearn.Models.Distributions.Multivariate.Internal.Container +import HLearn.Models.Distributions.Multivariate.Internal.Ignore +import HLearn.Models.Distributions.Multivariate.Internal.Marginalization +import HLearn.Models.Distributions.Multivariate.Internal.Unital +import HLearn.Models.Distributions.Multivariate.Internal.TypeLens +import HLearn.Models.Distributions.Multivariate.MultiNormal + +------------------------------------------------------------------------------- +-- Multivariate + +-- | this is the main type for specifying multivariate distributions +newtype Multivariate (dp:: *) (xs :: [[* -> * -> *]]) prob = Multivariate (MultivariateTF (Concat xs) prob) + +type family MultivariateTF (xs::[* -> * -> *]) prob +type instance MultivariateTF '[] prob = Unital prob +type instance MultivariateTF ((Container univariate sample) ': xs) prob = + Container univariate sample (MultivariateTF xs prob) prob +type instance MultivariateTF ((MultiContainer dist sample) ': xs) prob = + MultiContainer dist sample (MultivariateTF xs prob) prob +type instance MultivariateTF ((CatContainer label) ': xs) prob = + CatContainer label (MultivariateTF xs prob) prob +type instance MultivariateTF ((Ignore' label) ': xs) prob = + Ignore' label (MultivariateTF xs prob) prob + +deriving instance (Read (MultivariateTF (Concat xs) prob)) => Read (Multivariate dp xs prob) +deriving instance (Show (MultivariateTF (Concat xs) prob)) => Show (Multivariate dp xs prob) +deriving instance (Eq (MultivariateTF (Concat xs) prob)) => Eq (Multivariate dp xs prob) +deriving instance (Ord (MultivariateTF (Concat xs) prob)) => Ord (Multivariate dp xs prob) +deriving instance (Monoid (MultivariateTF (Concat xs) prob)) => Monoid (Multivariate dp xs prob) +deriving instance (Group (MultivariateTF (Concat xs) prob)) => Group (Multivariate dp xs prob) +deriving instance (NFData (MultivariateTF (Concat xs) prob)) => NFData (Multivariate dp xs prob) + +instance + ( HomTrainer (MultivariateTF (Concat xs) prob) + , Trainable dp + , GetHList dp ~ Datapoint (MultivariateTF (Concat xs) prob) + ) => HomTrainer (Multivariate dp xs prob) + where + type Datapoint (Multivariate dp xs prob) = dp + train1dp dp = Multivariate $ train1dp $ getHList dp + +instance Probabilistic (Multivariate dp xs prob) where + type Probability (Multivariate dp xs prob) = prob + +instance + ( PDF (MultivariateTF (Concat xs) prob) + , Probability (MultivariateTF (Concat xs) prob) ~ prob + , Datapoint (MultivariateTF (Concat xs) prob) ~ GetHList dp + , Trainable dp + , HomTrainer (Multivariate dp xs prob) + ) => PDF (Multivariate dp xs prob) + where + pdf (Multivariate dist) dp = pdf dist (getHList dp) + +instance + ( Marginalize' (Nat1Box n) (MultivariateTF (Concat xs) prob) + , MarginalizeOut' (Nat1Box n) (MultivariateTF (Concat xs) prob) + ~ MultivariateTF (Concat (Replace2D n xs (Ignore' (Index (HList2TypeList (GetHList dp)) n)))) prob + ) => Marginalize' (Nat1Box n) (Multivariate dp xs prob) + where + + type Margin' (Nat1Box n) (Multivariate dp xs prob) = Margin' (Nat1Box n) (MultivariateTF (Concat xs) prob) + getMargin' n (Multivariate dist) = getMargin' n dist + + type MarginalizeOut' (Nat1Box n) (Multivariate dp xs prob) = + Multivariate dp (Replace2D n xs (Ignore' (Index (HList2TypeList (GetHList dp)) n))) prob + marginalizeOut' n (Multivariate dist) = Multivariate $ marginalizeOut' n dist + + condition' n (Multivariate dist) dp = Multivariate $ condition' n dist dp + +type family HList2TypeList hlist :: [a] +type instance HList2TypeList (HList xs) = xs + +type family Index (xs::[a]) (i::Nat1) :: a +type instance Index (x ': xs) Zero = x +type instance Index (x ': xs) (Succ i) = Index xs i + +type family Replace2D (n :: Nat1) (xs :: [ [ a ] ]) (newval :: a) :: [ [ a ] ] +type instance Replace2D Zero ((x ': xs) ': ys) newval = (newval ': xs) ': ys +type instance Replace2D (Succ n) ((x ': xs) ': ys) newval = AppendFront x (Replace2D n (xs ': ys) newval) +type instance Replace2D n ('[] ': ys) newval = '[] ': (Replace2D n ys newval) + +type family AppendFront (x :: a) (xs :: [[a]]) :: [[a]] +type instance AppendFront x (xs ': ys) = (x ': xs) ': ys + +data Boxer xs = Boxer + +------------------------------------------------------------------------------- +-- Type functions + +-- type Multivariate (xs::[[* -> * -> *]]) prob = MultivariateTF (Concat xs) prob + +type family MultiCategorical (xs :: [*]) :: [* -> * -> *] +type instance MultiCategorical '[] = ('[]) +type instance MultiCategorical (x ': xs) = (CatContainer x) ': (MultiCategorical xs) + +-- type Dependent dist (xs :: [*]) = '[ MultiContainer (dist xs) xs ] +type family Dependent (dist::a) (xs :: [*]) :: [* -> * -> *] +type instance Dependent dist xs = '[ MultiContainer (dist xs) xs ] + +type family Independent (dist :: a) (sampleL :: [*]) :: [* -> * -> *] +type instance Independent dist '[] = '[] +type instance Independent (dist :: * -> *) (x ': xs) = (Container dist x) ': (Independent dist xs) +type instance Independent (dist :: * -> * -> *) (x ': xs) = (Container (dist x) x) ': (Independent dist xs) + diff --git a/src/HLearn/Models/Distributions/Multivariate/Internal/CatContainer.hs b/src/HLearn/Models/Distributions/Multivariate/Internal/CatContainer.hs new file mode 100644 --- /dev/null +++ b/src/HLearn/Models/Distributions/Multivariate/Internal/CatContainer.hs @@ -0,0 +1,235 @@ +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE FunctionalDependencies #-} +{-# LANGUAGE PolyKinds #-} +{-# LANGUAGE StandaloneDeriving #-} + +-- | The categorical distribution is used for discrete data. It is also sometimes called the discrete distribution or the multinomial distribution. For more, see the wikipedia entry: +module HLearn.Models.Distributions.Multivariate.Internal.CatContainer +{- ( + -- * Data types + CatContainer (CatContainer) + , CatContainerParams (..) + + -- * Helper functions + , dist2list + , mostLikely + )-} + where + +import Control.DeepSeq +import Control.Monad.Random +import Data.List +import Data.List.Extras +import Debug.Trace + +import qualified Data.Map.Strict as Map +import qualified Data.Foldable as F + +import HLearn.Algebra +import HLearn.Models.Distributions.Common +import HLearn.Models.Distributions.Multivariate.Internal.Ignore +import HLearn.Models.Distributions.Multivariate.Internal.Unital +import HLearn.Models.Distributions.Multivariate.Internal.Marginalization +import HLearn.Models.Distributions.Univariate.Categorical + +------------------------------------------------------------------------------- +-- data types + +data CatContainer label basedist prob = CatContainer + { pdfmap :: !(Map.Map label basedist) + , probmap :: !(Map.Map label prob) + , catnumdp :: prob + } + deriving (Show,Read,Eq,Ord) + +-- instance (Show basedist, Show label, Show prob) => Show (CatContainer label basedist prob) where +-- show dist = "CatContainer " +-- -- ++"{ "++"params="++show (params dist) +-- ++"{ "++"pdfmap="++show (pdfmap dist) +-- ++", catnumdp="++show (catnumdp dist) +-- ++"}" + +instance (NFData label, NFData prob, NFData basedist) => + NFData (CatContainer label basedist prob) + where + rnf d = rnf $ pdfmap d + +-- type CatContainer label basedist prob = RegSG2Group (CatContainer label basedist prob) + +------------------------------------------------------------------------------- +-- Algebra + +instance (Ord label, Num prob, Monoid basedist) => Abelian (CatContainer label basedist prob) +instance (Ord label, Num prob, Monoid basedist) => Monoid (CatContainer label basedist prob) where + mempty = CatContainer mempty mempty 0 + d1 `mappend` d2 = CatContainer + { pdfmap = Map.unionWith (<>) (pdfmap d1) (pdfmap d2) + , probmap = Map.unionWith (+) (probmap d1) (probmap d2) + , catnumdp = (catnumdp d1)+(catnumdp d2) + } + +instance (Ord label, Num prob, Group basedist) => Group (CatContainer label basedist prob) where + inverse d1 = CatContainer + { pdfmap = Map.map (inverse) (pdfmap d1) + , probmap = Map.map negate (probmap d1) + , catnumdp = -catnumdp d1 + } + +instance (Num prob) => HasRing (CatContainer label basedist prob) where + type Ring (CatContainer label basedist prob) = prob + +instance + ( Ord label + , Num prob + , Module basedist + , Ring basedist ~ Ring (CatContainer label basedist prob) + ) => Module (CatContainer label basedist prob) + where + r .* d = CatContainer + { pdfmap = Map.map (r.*) (pdfmap d) + , probmap = Map.map (r*) (probmap d) + , catnumdp = r * catnumdp d + } +-- -- instance (Ord label, Num prob) => LeftModule prob (CatContainer label prob) +-- instance (Ord label, Num prob) => LeftOperator prob (CatContainer label prob) where +-- p .* (CatContainer pdf) = CatContainer $ Map.map (*p) pdf +-- +-- -- instance (Ord label, Num prob) => RightModule prob (CatContainer label prob) +-- instance (Ord label, Num prob) => RightOperator prob (CatContainer label prob) where +-- (*.) = flip (.*) + +------------------------------------------------------------------------------- +-- Training + +instance + ( Ord label + , Num prob + , HomTrainer basedist + , Datapoint basedist ~ HList ys + ) => HomTrainer (CatContainer label basedist prob) + where + type Datapoint (CatContainer label basedist prob) = label `HCons` (Datapoint basedist) + + train1dp (dp:::basedp) = CatContainer + { pdfmap = Map.singleton dp $ train1dp basedp + , probmap = Map.singleton dp 1 + , catnumdp = 1 + } + +instance (Num prob) => NumDP (CatContainer label basedist prob) where + numdp dist = catnumdp dist + +------------------------------------------------------------------------------- +-- Distribution + +instance Probabilistic (CatContainer label basedist prob) where + type Probability (CatContainer label basedist prob) = prob + +instance + ( Ord prob, Fractional prob, Show prob, Probability basedist ~ prob + , Ord label + , PDF basedist + , Datapoint basedist ~ HList ys + , Show (Datapoint basedist) + , Show label + ) => PDF (CatContainer label basedist prob) + where + + {-# INLINE pdf #-} + pdf dist (label:::basedp) = val*weight/(catnumdp dist) + where + weight = case Map.lookup label (probmap dist) of + Nothing -> 0 + Just x -> x + val = case Map.lookup label (pdfmap dist) of + Nothing -> trace ("Warning.CatContainer: label "++show label++" not found in training data: "++show (Map.keys $ pdfmap dist)) $ 0 + Just x -> pdf x basedp + +--------------------------------------- + +instance + ( NumDP basedist + , Ring basedist ~ prob + , Monoid basedist + , HCons label (Datapoint basedist) ~ HList (label ': ts) + , Ord label + ) => Marginalize' (Nat1Box Zero) (CatContainer label basedist prob) + where + + type Margin' (Nat1Box Zero) (CatContainer label basedist prob) = (Categorical label prob) + getMargin' _ dist = Categorical $ probmap dist --Map.map numdp (pdfmap dist) + + type MarginalizeOut' (Nat1Box Zero) (CatContainer label basedist prob) = Ignore' label basedist prob + marginalizeOut' _ dist = Ignore' $ reduce $ Map.elems (pdfmap dist) + + condition' _ dist dp = Ignore' $ + case Map.lookup dp (pdfmap dist) of + Nothing -> error "CatContainer.condition: Nothing" + Just basedist -> basedist + +{- conditionAllButOne _ dist (dp:::dpL) = Ignore' $ + case Map.lookup dp (pdfmap dist) of + Nothing -> error "CatContainer.condition: Nothing" + Just basedist -> basedist-} + + +instance + ( Marginalize' (Nat1Box n) basedist + , Monoid basedist + , PDF (Margin' (Nat1Box n) basedist) + , prob ~ Probability (Margin' (Nat1Box n) basedist) + , prob ~ Ring basedist + , Module basedist + , Ord label + , Num prob + ) => Marginalize' (Nat1Box (Succ n)) (CatContainer label basedist prob) + where + + type Margin' (Nat1Box (Succ n)) (CatContainer label basedist prob) = Margin' (Nat1Box n) basedist + getMargin' _ dist = getMargin' (undefined :: Nat1Box n) $ reduce $ + zipWith (.*) + (Map.elems $ probmap dist) + (Map.elems $ pdfmap dist) + + type MarginalizeOut' (Nat1Box (Succ n)) (CatContainer label basedist prob) = + CatContainer label (MarginalizeOut' (Nat1Box n) basedist) prob + marginalizeOut' _ dist = dist { pdfmap = fmap (marginalizeOut' (undefined :: Nat1Box n)) $ pdfmap dist } + + condition' _ dist dp = dist + { probmap = Map.unionWith (*) (probmap dist) (conditionmap) + , pdfmap = fmap (flip (condition' (undefined :: Nat1Box n)) dp) $ pdfmap dist + } + where + conditionmap = fmap (\dist -> pdf (getMargin' (undefined :: Nat1Box n) dist) dp) $ pdfmap dist +-- conditionAllButOne _ dist (dp:::dpL) = dist { pdfmap = fmap (flip (condition (undefined :: Nat1Box n)) dpL) $ pdfmap dist } + +{-marginalizeRight :: + ( NumDP basedist prob + ) => CatContainer label basedist prob -> CatContainer label (Unital prob) prob +marginalizeRight dist = CatContainer + { pdfmap = Map.map (Unital . numdp) (pdfmap dist) + , probmap = error "probmap" + , catnumdp = catnumdp dist + }-} +-- marginalizeRight (SGJust dist) = Map.foldr mappend mempty (pdfmap dist) + + +------------------------------------------------------------------------------- +-- test + +ds= [ "test":::'g':::"foo":::HNil + , "test":::'f':::"fok":::HNil + , "toot":::'f':::"foo":::HNil + ] + +test = train ds :: CatContainer String (CatContainer Char (CatContainer String (Unital Double) Double) Double) Double diff --git a/src/HLearn/Models/Distributions/Multivariate/Internal/Container.hs b/src/HLearn/Models/Distributions/Multivariate/Internal/Container.hs new file mode 100644 --- /dev/null +++ b/src/HLearn/Models/Distributions/Multivariate/Internal/Container.hs @@ -0,0 +1,228 @@ +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE FunctionalDependencies #-} +{-# LANGUAGE PolyKinds #-} +{-# LANGUAGE StandaloneDeriving #-} + +-- | +module HLearn.Models.Distributions.Multivariate.Internal.Container + ( Container + , MultiContainer + ) + where + +import Control.DeepSeq +import Debug.Trace +import GHC.TypeLits +import HLearn.Algebra +import HLearn.Models.Distributions.Common +import HLearn.Models.Distributions.Multivariate.Internal.Ignore +import HLearn.Models.Distributions.Multivariate.Internal.Marginalization + +------------------------------------------------------------------------------- +-- data types + +data Container dist sample basedist (prob:: * ) = Container + { dist :: dist prob + , basedist :: basedist + } + deriving (Read,Show,Eq,Ord) + +instance (NFData (dist prob), NFData basedist) => NFData (Container dist sample basedist prob) where + rnf c = deepseq (dist c) $ rnf (basedist c) + +newtype MultiContainer dist sample basedist prob = MultiContainer (Container dist sample basedist prob) + deriving (Read,Show,Eq,Ord,Monoid,Abelian,Group,NFData) + +------------------------------------------------------------------------------- +-- Algebra + +instance (Abelian (dist prob), Abelian basedist) => Abelian (Container dist sample basedist prob) +instance + ( Monoid (dist prob) + , Monoid basedist + ) => Monoid (Container dist sample basedist prob) + where + mempty = Container mempty mempty + c1 `mappend` c2 = Container + { dist = dist c1 <> dist c2 + , basedist = basedist c1 <> basedist c2 + } + +instance + ( Group (dist prob) + , Group basedist + ) => Group (Container dist sample basedist prob) + where + inverse c = Container + { dist = inverse $ dist c + , basedist = inverse $ basedist c + } + +instance + ( HasRing (dist prob) + , HasRing basedist + , Ring (dist prob) ~ Ring basedist + ) => HasRing (Container dist sample basedist prob) + where + type Ring (Container dist sample basedist prob) = Ring (dist prob) + + +instance + ( HasRing (dist prob) + , HasRing basedist + , Ring (dist prob) ~ Ring basedist + ) => HasRing (MultiContainer dist sample basedist prob) + where + type Ring (MultiContainer dist sample basedist prob) = Ring (dist prob) + +instance + ( Module (dist prob) + , Module basedist + , Ring (dist prob) ~ Ring basedist + ) => Module (Container dist sample basedist prob) + where + r .* c = Container + { dist = r .* (dist c) + , basedist = r .* (basedist c) + } + +deriving instance + ( Module (dist prob) + , Module basedist + , Ring (dist prob) ~ Ring basedist + ) => Module (MultiContainer dist sample basedist prob) + + +------------------------------------------------------------------------------- +-- Training + +instance + ( HomTrainer (dist prob) + , HomTrainer basedist + , Datapoint basedist ~ HList ys + ) => HomTrainer (Container dist sample basedist prob) + where + type Datapoint (Container dist sample basedist prob) = + (Datapoint (dist prob)) `HCons` (Datapoint basedist) + + train1dp (dp:::basedp) = Container + { dist = train1dp dp + , basedist = train1dp basedp + } + +instance (NumDP (dist prob), HasRing basedist, Ring basedist ~ Ring (dist prob)) => NumDP (Container dist sample basedist prob) where + numdp container = numdp $ dist container + +--------------------------------------- + +instance + ( HomTrainer (dist prob) + , HomTrainer basedist + , Datapoint (dist prob) ~ HList zs + , Datapoint basedist ~ HList ys + , HTake1 (Nat1Box (Length1 zs)) (HList (zs++ys)) (HList zs) + , HDrop1 (Nat1Box (Length1 zs)) (HList (zs++ys)) (HList ys) + ) => HomTrainer (MultiContainer dist sample basedist prob) + where + type Datapoint (MultiContainer dist sample basedist prob) = + (Datapoint (dist prob)) `HAppend` (Datapoint basedist) + + train1dp dpL = MultiContainer $ Container + { dist = train1dp $ htake1 (Nat1Box :: Nat1Box (Length1 zs)) dpL + , basedist = train1dp $ hdrop1 (Nat1Box :: Nat1Box (Length1 zs)) dpL + } + +instance (NumDP (dist prob), HasRing basedist, Ring basedist ~ Ring (dist prob)) => NumDP (MultiContainer dist sample basedist prob) where + numdp (MultiContainer container) = numdp $ dist container + +------------------------------------------------------------------------------- +-- Distribution + +instance Probabilistic (Container dist sample basedist prob) where + type Probability (Container dist sample basedist prob) = prob + +instance + ( PDF (dist prob) + , PDF basedist + , Probability (dist prob) ~ prob + , Probability basedist ~ prob + , Probabilistic (Container dist sample basedist prob) + , Datapoint basedist ~ HList ys + , Datapoint (dist prob) ~ y + , Datapoint (Container dist sample basedist prob) ~ HList (y ': ys) + , Num prob + ) => PDF (Container dist sample basedist prob) + where + pdf container (dp:::basedp) = pdf1*pdf2 + where + pdf1 = pdf (dist container) dp + pdf2 = pdf (basedist container) basedp + +instance Marginalize' (Nat1Box Zero) (Container dist (sample :: *) basedist prob) where + type Margin' (Nat1Box Zero) (Container dist sample basedist prob) = dist prob + getMargin' _ container = dist container + + type MarginalizeOut' (Nat1Box Zero) (Container dist sample basedist prob) = Ignore' sample basedist prob + marginalizeOut' _ container = Ignore' $ basedist container + + condition' _ container dp = Ignore' $ basedist container --error "Container.Marginalize.condition: undefined" + +instance + ( Marginalize' (Nat1Box n) basedist + ) => Marginalize' (Nat1Box (Succ n)) (Container dist sample basedist prob) + where + type Margin' (Nat1Box (Succ n)) (Container dist sample basedist prob) = Margin' (Nat1Box n) basedist + getMargin' _ container = getMargin' (undefined :: Nat1Box n) $ basedist container + + type MarginalizeOut' (Nat1Box (Succ n)) (Container dist sample basedist prob) = + Container dist sample (MarginalizeOut' (Nat1Box n) basedist) prob + marginalizeOut' _ container = Container + { dist = dist container + , basedist = marginalizeOut' (undefined :: Nat1Box n) $ basedist container + } + + condition' _ container dp = Container + { dist = dist container + , basedist = condition' (undefined :: Nat1Box n) (basedist container) dp + } + +{-instance Marginalize (Nat1Box Zero) (Container dist sample basedist prob) (dist prob) where + getMargin _ container = dist container + +instance + ( Marginalize (Nat1Box n) basedist margin + ) => Marginalize (Nat1Box (Succ n)) (Container dist sample basedist prob) margin + where + getMargin _ container = getMargin (undefined :: Nat1Box n) $ basedist container +-} +--------------------------------------- + +instance Probabilistic (MultiContainer dist sample basedist prob) where + type Probability (MultiContainer dist sample basedist prob) = prob + +instance + ( PDF (dist prob) + , PDF basedist + , prob ~ Probability (dist prob) + , prob ~ Probability basedist + , Num prob + , Datapoint (dist prob) ~ HList dpL + , Datapoint basedist ~ HList basedpL + , HTake1 (Nat1Box (Length1 dpL)) (HList (dpL ++ basedpL)) (HList dpL) + , HDrop1 (Nat1Box (Length1 dpL)) (HList (dpL ++ basedpL)) (HList basedpL) + ) => PDF (MultiContainer dist sample basedist prob) + where + pdf (MultiContainer container) dp = (pdf (dist container) dp1)*(pdf (basedist container) dp2) + where + dp1 = htake1 (Nat1Box :: Nat1Box (Length1 dpL)) dp + dp2 = hdrop1 (Nat1Box :: Nat1Box (Length1 dpL)) dp diff --git a/src/HLearn/Models/Distributions/Multivariate/Internal/Ignore.hs b/src/HLearn/Models/Distributions/Multivariate/Internal/Ignore.hs new file mode 100644 --- /dev/null +++ b/src/HLearn/Models/Distributions/Multivariate/Internal/Ignore.hs @@ -0,0 +1,103 @@ +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE FunctionalDependencies #-} +{-# LANGUAGE PolyKinds #-} +{-# LANGUAGE StandaloneDeriving #-} + +-- | Used for ignoring data +module HLearn.Models.Distributions.Multivariate.Internal.Ignore + ( Ignore + , Ignore' (Ignore') + ) where + +import Control.DeepSeq + +import HLearn.Algebra +import HLearn.Models.Distributions.Common +import HLearn.Models.Distributions.Multivariate.Internal.Unital +import HLearn.Models.Distributions.Multivariate.Internal.Marginalization + +------------------------------------------------------------------------------- +-- data types + +newtype Ignore' (label:: *) (basedist:: *) (prob :: *) = Ignore' { basedist :: basedist } + deriving (Show,Read,Eq,Ord,NFData) + +type family Ignore (xs :: [*]) :: [* -> * -> *] +type instance Ignore '[] = '[] +type instance Ignore (x ': xs) = (Ignore' x) ': (Ignore xs) + +------------------------------------------------------------------------------- +-- Algebra + +instance (Abelian basedist) => Abelian (Ignore' label basedist prob) +instance (Monoid basedist) => Monoid (Ignore' label basedist prob) where + mempty = Ignore' mempty + mappend d1 d2 = Ignore' $ mappend (basedist d1) (basedist d2) + +instance (Group basedist) => Group (Ignore' label basedist prob) where + inverse d = Ignore' $ inverse (basedist d) + +instance (HasRing basedist) => HasRing (Ignore' label basedist prob) where + type Ring (Ignore' label basedist prob) = Ring basedist + +instance (Module basedist) => Module (Ignore' label basedist prob) where + r .* d = Ignore' $ r .* basedist d + +------------------------------------------------------------------------------- +-- Training + +instance + ( HomTrainer basedist + , Datapoint basedist ~ HList ys + ) => HomTrainer (Ignore' label basedist prob) + where + type Datapoint (Ignore' label basedist prob) = label `HCons` (Datapoint basedist) + + train1dp (dp:::basedp) = Ignore' $ train1dp basedp + +instance (NumDP basedist) => NumDP (Ignore' label basedist prob) where + numdp (Ignore' basedist) = numdp basedist + +------------------------------------------------------------------------------- +-- Distribution + +instance Probabilistic (Ignore' label basedist prob) where + type Probability (Ignore' label basedist prob) = prob + +instance + ( Probability basedist ~ prob + , HomTrainer (Ignore' label basedist prob) + , Datapoint (Ignore' label basedist prob) ~ HList dpL + , Datapoint basedist ~ HList basedpL + , PDF basedist + ) => PDF (Ignore' label basedist prob) + where + + {-# INLINE pdf #-} + pdf dist (label:::basedp) = pdf (basedist dist) basedp + +-- instance Marginalize (Nat1Box Zero) (Ignore' label basedist prob) (Unital prob) where +-- getMargin _ dist = Categorical $ Map.map numdp (pdfmap dist) + +instance + ( Marginalize' (Nat1Box n) basedist + ) => Marginalize' (Nat1Box (Succ n)) (Ignore' label basedist prob) + where + type Margin' (Nat1Box (Succ n)) (Ignore' label basedist prob) = Margin' (Nat1Box n) basedist + getMargin' _ dist = getMargin' (undefined :: Nat1Box n) $ basedist dist + + type MarginalizeOut' (Nat1Box (Succ n)) (Ignore' label basedist prob) = + Ignore' label (MarginalizeOut' (Nat1Box n) basedist) prob + marginalizeOut' _ dist = Ignore' $ marginalizeOut' (undefined :: Nat1Box n) $ basedist dist + + condition' _ dist dp = Ignore' $ condition' (undefined :: Nat1Box n) (basedist dist) dp diff --git a/src/HLearn/Models/Distributions/Multivariate/Internal/Marginalization.hs b/src/HLearn/Models/Distributions/Multivariate/Internal/Marginalization.hs new file mode 100644 --- /dev/null +++ b/src/HLearn/Models/Distributions/Multivariate/Internal/Marginalization.hs @@ -0,0 +1,66 @@ +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE FunctionalDependencies #-} +{-# LANGUAGE PolyKinds #-} +{-# LANGUAGE StandaloneDeriving #-} + +-- | Marginalization is a tricky procedure involving a lot of type hackery. All user-facing interfaces use the "Marginalize" class, and all internal interfaces use the "Marginalize'" class. Essentially, "TypeLens" class converts the "Type lens" for our user's data type into a type level number. The "Marginalize'" class then takes this type number and uses it as an index into an appropriate type list that represents the data type. +-- +-- The TemplateHaskell module has more info. + +module HLearn.Models.Distributions.Multivariate.Internal.Marginalization + where + +import GHC.TypeLits +import HLearn.Algebra +import HLearn.Models.Distributions.Common +import HLearn.Models.Distributions.Multivariate.Internal.TypeLens + +------------------------------------------------------------------------------- +-- external type classes + +class (Marginalize' (TypeLensIndex index) dist, TypeLens index) => Marginalize index dist where + + type Margin index dist + getMargin :: index -> dist -> Margin index dist + + type MarginalizeOut index dist + marginalizeOut :: index -> dist -> MarginalizeOut index dist + + condition :: index -> Datapoint (Margin' (TypeLensIndex index) dist) -> dist -> MarginalizeOut' (TypeLensIndex index) dist + +instance + ( Marginalize' (TypeLensIndex index) dist + , TypeLens index + ) => Marginalize index dist + where + + type Margin index dist = Margin' (TypeLensIndex index) dist + getMargin _ dist = getMargin' (undefined :: (TypeLensIndex index)) dist + + type MarginalizeOut index dist = MarginalizeOut' (TypeLensIndex index) dist + marginalizeOut _ dist = marginalizeOut' (undefined :: (TypeLensIndex index)) dist + condition _ dp dist = condition' (undefined :: (TypeLensIndex index)) dist dp + +------------------------------------------------------------------------------- +-- internal type classes + +class Marginalize' index dist where + type Margin' index dist + getMargin' :: index -> dist -> Margin' index dist + + type MarginalizeOut' index dist + marginalizeOut' :: index -> dist -> MarginalizeOut' index dist + + condition' :: index -> dist -> Datapoint (Margin' index dist) -> MarginalizeOut' index dist + +-- conditionAllButOne :: index -> dist -> Datapoint dist -> MarginalizeOut index dist diff --git a/src/HLearn/Models/Distributions/Multivariate/Internal/TypeLens.hs b/src/HLearn/Models/Distributions/Multivariate/Internal/TypeLens.hs new file mode 100644 --- /dev/null +++ b/src/HLearn/Models/Distributions/Multivariate/Internal/TypeLens.hs @@ -0,0 +1,196 @@ +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE FlexibleContexts #-} + +-- | This module provides convenient TemplateHaskell functions for making type lens suitable for use with multivariate distributions. +-- +-- Given a data type that looks like: +-- +-- >data Character = Character +-- > { _name :: String +-- > , _species :: String +-- > , _job :: Job +-- > , _isGood :: Maybe Bool +-- > , _age :: Double -- in years +-- > , _height :: Double -- in feet +-- > , _weight :: Double -- in pounds +-- > } +-- > deriving (Read,Show,Eq,Ord) +-- > +-- >data Job = Manager | Crew | Henchman | Other +-- > deriving (Read,Show,Eq,Ord) +-- +-- +-- when we run the command: +-- +-- >makeTypeLenses ''Character +-- +-- We generate the following type lenses automatically: +-- +-- >data TH_name = TH_name +-- >data TH_species = TH_species +-- >data TH_job = TH_job +-- >data TH_isGood = TH_isGood +-- >data TH_age = TH_age +-- >data TH_height = TH_height +-- >data TH_weight = TH_weight +-- > +-- >instance TypeLens TH_name where +-- > type instance TypeLensIndex TH_name = Nat1Box Zero +-- >instance TypeLens TH_species where +-- > type instance TypeLensIndex TH_species = Nat1Box (Succ Zero) +-- >instance TypeLens TH_job where +-- > type instance TypeLensIndex TH_job = Nat1Box (Succ (Succ Zero)) +-- >instance TypeLens TH_isGood where +-- > type instance TypeLensIndex TH_isGood = Nat1Box (Succ (Succ (Succ Zero))) +-- >instance TypeLens TH_age where +-- > type instance TypeLensIndex TH_age = Nat1Box (Succ (Succ (Succ (Succ Zero)))) +-- >instance TypeLens TH_height where +-- > type instance TypeLensIndex TH_height = Nat1Box (Succ (Succ (Succ (Succ (Succ Zero))))) +-- >instance TypeLens TH_weight where +-- > type instance TypeLensIndex TH_weight = Nat1Box (Succ (Succ (Succ (Succ (Succ (Succ Zero)))))) +-- > +-- >instance Trainable Character where +-- > type instance GetHList Character = HList '[String,String,Job,Maybe Bool, Double,Double,Double] +-- > getHList var = name var:::species var:::job var:::isGood var:::age var:::height var:::weight var:::HNil +-- > +-- >instance MultivariateLabels Character where +-- > getLabels dist = ["TH_name","TH_species","TH_job","TH_isGood","TH_age","TH_height","TH_weight"] +-- +-- +-- + +module HLearn.Models.Distributions.Multivariate.Internal.TypeLens + ( + -- * Lens + Trainable (..) + , TypeLens (..) + -- * TemplateHaskell + , makeTypeLenses + , nameTransform + ) + where + +import HLearn.Algebra +import Language.Haskell.TH +import Language.Haskell.TH.Syntax + + +------------------------------------------------------------------------------- +-- Trainable + +-- | The Trainable class allows us to convert data types into an isomorphic "HList". All of our multivariate distributions work on "HList"s, so they work on all instances of "Trainable" as well. +class Trainable t where + type GetHList t + getHList :: t -> GetHList t + +instance Trainable (HList '[]) where + type GetHList (HList '[]) = HList '[] + getHList t = t + +instance (Trainable (HList xs)) => Trainable (HList (x ': xs)) where + type GetHList (HList (x ': xs)) = HList (x ': xs) + getHList t = t + +-- | This specifies a type level natural number (i.e. "Nat1") that indexes at the right location into our HList +class TypeLens i where + type TypeLensIndex i + +-- | given the name of one of our records, transform it into the name of our type lens +nameTransform :: String -> String +nameTransform str = "TH"++str + +-- | constructs the type lens +makeTypeLenses :: Name -> Q [Dec] +makeTypeLenses name = do + datatypes <- makeDatatypes name + indexNames <- makeIndexNames name + trainableInstance <- makeTrainable name + multivariateLabels <- makeMultivariateLabels name + return $ datatypes ++ indexNames ++ trainableInstance ++ multivariateLabels + +makeDatatypes :: Name -> Q [Dec] +makeDatatypes name = fmap (map makeEmptyData) $ extractContructorNames name + where + makeEmptyData str = DataD [] (mkName $ nameTransform str) [] [NormalC (mkName $ nameTransform str) []] [] + +makeIndexNames :: Name -> Q [Dec] +makeIndexNames name = fmap (map makeIndexName . zip [0..]) $ extractContructorNames name + where + makeIndexName (i,str) = InstanceD [] (AppT (ConT $ mkName "TypeLens") (ConT $ mkName $ nameTransform str)) + [ TySynInstD (mkName "TypeLensIndex") [ConT $ mkName $ nameTransform str] (AppT (ConT $ mkName "Nat1Box") (typeNat i)) + ] + where + typeNat 0 = ConT $ mkName "Zero" + typeNat n = AppT (ConT $ mkName "Succ") $ typeNat (n-1) + +makeTrainable :: Name -> Q [Dec] +makeTrainable name = do + hlistType <- extractHListType name + hlistExp <- extractHListExp (mkName "var") name + return $ [InstanceD [] (AppT (ConT (mkName "Trainable")) (ConT name)) + [ TySynInstD (mkName "GetHList") [ConT name] (AppT (ConT $ mkName "HList") hlistType) + , FunD (mkName "getHList") [Clause [VarP $ mkName "var"] (NormalB hlistExp) []] + ]] + +makeMultivariateLabels :: Name -> Q [Dec] +makeMultivariateLabels name = do + labelL <- extractContructorNames name + return $ [ InstanceD [] (AppT (ConT (mkName "MultivariateLabels")) (ConT name)) + [ FunD (mkName "getLabels") [Clause [VarP $ mkName "dist"] (NormalB $ go labelL ) []] + ]] + where + go [] = ConE $ mkName "[]" + go (x:xs) = AppE (AppE (ConE $ mkName ":") (LitE $ StringL (nameTransform x))) $ go xs + +------------------------------------------------------------------------------- +-- below taken from Data.Lens + +type ConstructorFieldInfo = (Name, Strict, Type) + +extractHListType :: Name -> Q Type +extractHListType name = do + typeL <- fmap (map getType) $ extractConstructorFields name + return $ go typeL + where + go [] = ConT $ mkName "[]" + go (x:xs) = AppT (AppT (ConT $ mkName ":") (x)) $ go xs + + getType (n,s,t) = t + +extractHListExp :: Name -> Name -> Q Exp +extractHListExp var name = do + typeL <- fmap (map getName) $ extractConstructorFields name + return $ go typeL + where + go [] = ConE $ mkName "HNil" + go (x:xs) = AppE (AppE (ConE $ mkName ":::") (AppE (VarE x) (VarE var))) $ go xs + + getName (n,s,t) = n + +extractContructorNames :: Name -> Q [String] +extractContructorNames datatype = fmap (map name) $ extractConstructorFields datatype + where + name (n,s,t) = nameBase n + +extractConstructorFields :: Name -> Q [ConstructorFieldInfo] +extractConstructorFields datatype = do + let datatypeStr = nameBase datatype + i <- reify datatype + return $ case i of + TyConI (DataD _ _ _ [RecC _ fs] _) -> fs + TyConI (NewtypeD _ _ _ (RecC _ fs) _) -> fs + TyConI (DataD _ _ _ [_] _) -> error $ "Can't derive Lens without record selectors: " ++ datatypeStr + TyConI NewtypeD{} -> error $ "Can't derive Lens without record selectors: " ++ datatypeStr + TyConI TySynD{} -> error $ "Can't derive Lens for type synonym: " ++ datatypeStr + TyConI DataD{} -> error $ "Can't derive Lens for tagged union: " ++ datatypeStr + _ -> error $ "Can't derive Lens for: " ++ datatypeStr ++ ", type name required." + +extractTypeInfo name = do + i <- reify name + return $ case i of + TyConI (DataD _ n ts _ _) -> (n, ts) + TyConI (NewtypeD _ n ts _ _) -> (n, ts) diff --git a/src/HLearn/Models/Distributions/Multivariate/Internal/Unital.hs b/src/HLearn/Models/Distributions/Multivariate/Internal/Unital.hs new file mode 100644 --- /dev/null +++ b/src/HLearn/Models/Distributions/Multivariate/Internal/Unital.hs @@ -0,0 +1,61 @@ +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} + +module HLearn.Models.Distributions.Multivariate.Internal.Unital + where + +import Control.DeepSeq + +import HLearn.Algebra +import HLearn.Models.Distributions.Common + +------------------------------------------------------------------------------- +-- data types + +newtype Unital prob = Unital prob + deriving (Read,Show,Eq,Ord,NFData) + +------------------------------------------------------------------------------- +-- algebra + +instance (Num prob) => Abelian (Unital prob) where +instance (Num prob) => Monoid (Unital prob) where + mempty = Unital 0 + (Unital p1) `mappend` (Unital p2) = Unital $ p1+p2 + +instance (Num prob) => Group (Unital prob) where + inverse (Unital p) = Unital (-p) + +instance (Num prob) => HasRing (Unital prob) where + type Ring (Unital prob) = prob + +instance (Num prob) => Module (Unital prob) where + r .* (Unital p) = Unital $ r*p + +------------------------------------------------------------------------------- +-- training + +instance (Num prob) => HomTrainer (Unital prob) where + type Datapoint (Unital prob) = HList '[] + train1dp HNil = Unital 1 + +instance (Num prob) => NumDP (Unital prob) where + numdp (Unital prob) = prob + +------------------------------------------------------------------------------- +-- distributions + +instance (Num prob) => Probabilistic (Unital prob) where + type Probability (Unital prob) = prob + +instance (Num prob) => PDF (Unital prob) where + pdf (Unital prob) HNil = 1 --prob + diff --git a/src/HLearn/Models/Distributions/Multivariate/MultiNormal.hs b/src/HLearn/Models/Distributions/Multivariate/MultiNormal.hs new file mode 100644 --- /dev/null +++ b/src/HLearn/Models/Distributions/Multivariate/MultiNormal.hs @@ -0,0 +1,217 @@ +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FunctionalDependencies #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE ConstraintKinds #-} +{-# LANGUAGE InstanceSigs #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE PolyKinds #-} +{-# LANGUAGE StandaloneDeriving #-} + +module HLearn.Models.Distributions.Multivariate.MultiNormal + ( MultiNormal (..) + ) + where + +-- import qualified Control.ConstraintKinds as CK +import qualified Data.Foldable as F +import qualified Data.Vector as V +import qualified Data.Vector.Unboxed as VU +import Data.Vector.Unboxed.Deriving + +import Control.DeepSeq +import GHC.ST +import GHC.TypeLits + +import Foreign.Storable +import Numeric.LinearAlgebra hiding ((<>)) +import qualified Numeric.LinearAlgebra as LA + +import HLearn.Algebra +import HLearn.Models.Distributions.Common + +------------------------------------------------------------------------------- +-- data types + +data MultiNormalVec (n::Nat) prob = MultiNormalVec + { q0 :: !prob + , q1 :: !(VU.Vector prob) + , q2 :: !(V.Vector (VU.Vector prob)) + } + deriving (Read,Show,Eq,Ord) + +instance NFData (MultiNormalVec n prob) where + rnf mn = seq mn () + +newtype MultiNormal (xs::[*]) prob = MultiNormal (MultiNormalVec (Length xs) prob) + deriving (Read,Show,Eq,Ord,NFData) + +deriving instance (Monoid (MultiNormalVec (Length xs) prob)) => Monoid (MultiNormal xs prob) +deriving instance (Abelian (MultiNormalVec (Length xs) prob)) => Abelian (MultiNormal xs prob) +deriving instance (Group (MultiNormalVec (Length xs) prob)) => Group (MultiNormal xs prob) +deriving instance (Module (MultiNormalVec (Length xs) prob)) => Module (MultiNormal xs prob) + +------------------------------------------------------------------------------- +-- algebra + +instance (Num prob, VU.Unbox prob, SingI n) => Abelian (MultiNormalVec n prob) +instance (Num prob, VU.Unbox prob, SingI n) => Monoid (MultiNormalVec n prob) where + mempty = MultiNormalVec + { q0 = 0 + , q1 = VU.replicate n 0 + , q2 = V.replicate n (VU.replicate n 0) + } + where + n = fromIntegral $ fromSing (sing :: Sing n) + mn1 `mappend` mn2 = MultiNormalVec + { q0 = (q0 mn1) + (q0 mn2) + , q1 = VU.zipWith (+) (q1 mn1) (q1 mn2) + , q2 = V.zipWith (VU.zipWith (+)) (q2 mn1) (q2 mn2) + } + +instance (Num prob, VU.Unbox prob, SingI n) => Group (MultiNormalVec n prob) where + inverse mn = MultiNormalVec + { q0 = negate $ q0 mn + , q1 = VU.map negate (q1 mn) + , q2 = V.map (VU.map negate) (q2 mn) + } + +instance (Num prob) => HasRing (MultiNormalVec n prob) where + type Ring (MultiNormalVec n prob) = prob + +instance (Num prob, VU.Unbox prob, SingI n) => Module (MultiNormalVec n prob) where + r .* mn = MultiNormalVec + { q0 = r * q0 mn + , q1 = VU.map (r*) (q1 mn) + , q2 = V.map (VU.map (r*)) (q2 mn) + } + +--------------------------------------- + +instance (Num prob) => HasRing (MultiNormal xs prob) where + type Ring (MultiNormal xs prob) = prob + +------------------------------------------------------------------------------- +-- training + +instance (SingI n, Num prob, VU.Unbox prob) => HomTrainer (MultiNormalVec n prob) where + type Datapoint (MultiNormalVec n prob) = (VU.Vector prob) + train1dp dp = MultiNormalVec + { q0 = 1 + , q1 = dp + , q2 = V.generate n (\i -> VU.generate n (\j -> (dp VU.! i)*(dp VU.! j))) + } + where + n = fromIntegral $ fromSing (sing :: Sing n) + +instance + ( SingI (Length xs) + , Num prob + , VU.Unbox prob + , HList2List (Datapoint (MultiNormal xs prob)) prob + ) => HomTrainer (MultiNormal xs prob) + where + type Datapoint (MultiNormal xs prob) = HList xs + train1dp dp = MultiNormal $ train1dp $ VU.fromList $ hlist2list dp + +instance (Num prob) => NumDP (MultiNormal xs prob) where + numdp (MultiNormal mn) = q0 mn + +------------------------------------------------------------------------------- +-- distribution + +class (Probabilistic dist) => Covariance dist where + covar :: dist -> Matrix (Probability dist) + +instance + ( VU.Unbox prob + , SingI k + , Num prob + ) => Probabilistic (MultiNormalVec k prob) + where + type Probability (MultiNormalVec k prob) = prob + +instance + ( VU.Unbox prob + , SingI k + , Fractional prob + , Enum prob + , Storable prob + ) => Covariance (MultiNormalVec k prob) + where + covar mn = (k> Probabilistic (MultiNormal dpL prob) + where + type Probability (MultiNormal dpL prob) = prob + +instance + ( HList2List (HList dpL) prob + , VU.Unbox prob + , Floating prob + , Field prob + , Enum prob + , SingI (FromNat1 (Length1 dpL)) +-- , Covariance (MultiNormal dpL prob) + , Storable prob + ) => PDF (MultiNormal dpL prob) + where + pdf (MultiNormal dist) dpL = 1/(sqrt $ (2*pi)^(k)*(det sigma))*(exp $ (-1/2)*(top) ) + where + top=minElement $ ((trans $ x `sub` mu) LA.<> (inv sigma) LA.<> (x `sub` mu)) + + k = fromIntegral $ fromSing (sing :: Sing (Length dpL)) :: Int + x = k><1 $ hlist2list dpL + n = q0 dist + sigma = covar dist + mu = k><1 $ map (/n) $ VU.toList (q1 $ dist) +-- covarM = (k> HomTrainer (Binomial sample prob) where + type Datapoint (Binomial sample prob) = sample + train1dp dp = Binomial $ train1dp dp + +------------------------------------------------------------------------------- +-- distribution + +instance (Num sample) => Probabilistic (Binomial sample prob) where + type Probability (Binomial sample prob) = prob + +instance (Floating prob) => PDF (Binomial Int Double) where + pdf (Binomial dist) dp = S.probability (binomial n p) dp + where + n = bin_n $ Binomial dist + p = bin_p $ Binomial dist + +bin_n :: Binomial Int Double -> Int +bin_n (Binomial dist) = round $ ((fromIntegral $ m1 dist :: Double) / (fromIntegral $ m0 dist)) / (bin_p $ Binomial dist) + +bin_p :: Binomial Int Double -> Double +bin_p (Binomial dist) = ((fromIntegral $ m1 dist) / (fromIntegral $ m0 dist)) + 1 - (fromIntegral $ m2 dist)/(fromIntegral $ m1 dist) + +instance + ( PDF (Binomial sample prob) +-- , PlottableDataPoint sample + , Show prob + , Show sample + , Ord sample + , Ord prob + , Num prob + , Integral sample + ) => PlottableDistribution (Binomial sample prob) +-- instance PlottableDistribution (Poisson Int Double) + where + + plotType _ = Points + + samplePoints dist = [min..max] + where + min = 0 + max = maximum [20,floor $ 3*mu] + mu = 5 + +-- instance (Fractional prob) => Mean (Binomial prob) where +-- mean (Binomial dist) = +-- +-- instance (Fractional prob) => Variance (Binomiral prob) where +-- variance dist = mean dist diff --git a/src/HLearn/Models/Distributions/Univariate/Categorical.hs b/src/HLearn/Models/Distributions/Univariate/Categorical.hs new file mode 100644 --- /dev/null +++ b/src/HLearn/Models/Distributions/Univariate/Categorical.hs @@ -0,0 +1,144 @@ +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE TypeFamilies #-} + + +-- | The categorical distribution is used for discrete data. It is also sometimes called the discrete distribution or the multinomial distribution. For more, see the wikipedia entry: +module HLearn.Models.Distributions.Univariate.Categorical + ( + -- * Data types + Categorical (Categorical) + + -- * Helper functions + , dist2list + , mostLikely + ) + where + +import Control.DeepSeq +import Control.Monad.Random +import Data.List +import Data.List.Extras +import Debug.Trace + +import qualified Data.Map.Strict as Map +import qualified Data.Foldable as F + +import HLearn.Algebra +import HLearn.Models.Distributions.Common +import HLearn.Models.Distributions.Visualization.Gnuplot + +------------------------------------------------------------------------------- +-- Categorical + +data Categorical sampletype prob = Categorical + { pdfmap :: !(Map.Map sampletype prob) + } + deriving (Show,Read,Eq,Ord) + +instance (NFData sampletype, NFData prob) => NFData (Categorical sampletype prob) where + rnf d = rnf $ pdfmap d + +------------------------------------------------------------------------------- +-- Algebra + +instance (Ord label, Num prob) => Abelian (Categorical label prob) +instance (Ord label, Num prob) => Monoid (Categorical label prob) where + mempty = Categorical Map.empty + mappend !d1 !d2 = Categorical $ res + where + res = Map.unionWith (+) (pdfmap d1) (pdfmap d2) + +instance (Ord label, Num prob) => Group (Categorical label prob) where + inverse d1 = d1 {pdfmap=Map.map (0-) (pdfmap d1)} + +instance (Num prob) => HasRing (Categorical label prob) where + type Ring (Categorical label prob) = prob +instance (Ord label, Num prob) => Module (Categorical label prob) where + p .* (Categorical pdf) = Categorical $ Map.map (*p) pdf + +------------------------------------------------------------------------------- +-- Training + +instance (Ord label, Num prob) => HomTrainer (Categorical label prob) where + type Datapoint (Categorical label prob) = label + train1dp dp = Categorical $ Map.singleton dp 1 + +instance (Num prob) => NumDP (Categorical label prob) where + numdp dist = F.foldl' (+) 0 $ pdfmap dist + +------------------------------------------------------------------------------- +-- Distribution + +instance Probabilistic (Categorical label prob) where + type Probability (Categorical label prob) = prob + +instance (Ord label, Ord prob, Fractional prob) => PDF (Categorical label prob) where + + {-# INLINE pdf #-} + pdf dist label = {-0.0001+-}(val/tot) + where + val = case Map.lookup label (pdfmap dist) of + Nothing -> 0 + Just x -> x + tot = F.foldl' (+) 0 $ pdfmap dist + +instance (Ord label, Ord prob, Fractional prob) => CDF (Categorical label prob) where + + {-# INLINE cdf #-} + cdf dist label = (Map.foldl' (+) 0 $ Map.filterWithKey (\k a -> k<=label) $ pdfmap dist) + / (Map.foldl' (+) 0 $ pdfmap dist) + + {-# INLINE cdfInverse #-} + cdfInverse dist prob = go cdfL + where + cdfL = sortBy (\(k1,p1) (k2,p2) -> compare p2 p1) $ map (\k -> (k,pdf dist k)) $ Map.keys $ pdfmap dist + go (x:[]) = fst $ last cdfL + go (x:xs) = if prob < snd x -- && prob > (snd $ head xs) + then fst x + else go xs +-- cdfInverse dist prob = argmax (cdf dist) $ Map.keys $ pdfmap dist + +-- {-# INLINE mean #-} +-- mean dist = fst $ argmax snd $ Map.toList $ pdfmap dist +-- +-- {-# INLINE drawSample #-} +-- drawSample dist = do +-- x <- getRandomR (0,1) +-- return $ cdfInverse dist (x::prob) + + +instance (Num prob, Ord prob, Ord label) => Mean (Categorical label prob) where + mean dist = fst $ argmax snd $ Map.toList $ pdfmap dist + +-- | Extracts the element in the distribution with the highest probability +mostLikely :: Ord prob => Categorical label prob -> label +mostLikely dist = fst $ argmax snd $ Map.toList $ pdfmap dist + +-- | Converts a distribution into a list of (sample,probability) pai +dist2list :: Categorical sampletype prob -> [(sampletype,prob)] +dist2list (Categorical pdfmap) = Map.toList pdfmap + + +instance + ( Ord label, Show label + , Ord prob, Show prob, Fractional prob + ) => PlottableDistribution (Categorical label prob) + where + samplePoints (Categorical dist) = Map.keys dist + plotType dist = Bar + +------------------------------------------------------------------------------- +-- Morphisms + +-- instance +-- ( Ord label +-- , Num prob +-- ) => Morphism (Categorical label prob) FreeModParams (FreeMod prob label) +-- where +-- Categorical pdf $> FreeModParams = FreeMod pdf +-- diff --git a/src/HLearn/Models/Distributions/Univariate/Exponential.hs b/src/HLearn/Models/Distributions/Univariate/Exponential.hs new file mode 100644 --- /dev/null +++ b/src/HLearn/Models/Distributions/Univariate/Exponential.hs @@ -0,0 +1,68 @@ +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE FunctionalDependencies #-} + +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE KindSignatures #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} + +-- | The method of moments can be used to estimate a number of commonly used distributions. This module is still under construction as I work out the best way to handle morphisms from the Moments3 type to types of other distributions. For more information, see the wikipedia entry: + +module HLearn.Models.Distributions.Univariate.Exponential + ( Exponential + ) + where + +import Control.DeepSeq +import GHC.TypeLits +import qualified Data.Vector.Unboxed as U +import Data.Vector.Unboxed.Deriving + +import HLearn.Algebra +import HLearn.Models.Distributions.Common +import HLearn.Models.Distributions.Univariate.Internal.Moments +import HLearn.Models.Distributions.Visualization.Gnuplot + +------------------------------------------------------------------------------- +-- Exponential + +newtype Exponential prob = Exponential (Moments3 prob) + deriving (Read,Show,Eq,Ord,Monoid,Group) + +instance (Num prob) => HomTrainer (Exponential prob) where + type Datapoint (Exponential prob) = prob + train1dp dp = Exponential $ train1dp dp + +instance (Num prob) => Probabilistic (Exponential prob) where + type Probability (Exponential prob) = prob + +instance (Floating prob) => PDF (Exponential prob) where + pdf dist dp = lambda*(exp $ (-1)*lambda*dp) + where + lambda = e_lambda dist + +e_lambda (Exponential dist) = (m0 dist)/(m1 dist) + +instance (Fractional prob) => Mean (Exponential prob) where + mean dist = 1/(e_lambda dist) + +instance (Fractional prob) => Variance (Exponential prob) where + variance dist = 1/(e_lambda dist)^^2 + +instance + ( Floating prob + , Enum prob + , Show prob + , Ord prob + ) => PlottableDistribution (Exponential prob) where + + plotType _ = Continuous + + samplePoints dist = samplesFromMinMax 0 $ (mean dist)*3 diff --git a/src/HLearn/Models/Distributions/Univariate/Geometric.hs b/src/HLearn/Models/Distributions/Univariate/Geometric.hs new file mode 100644 --- /dev/null +++ b/src/HLearn/Models/Distributions/Univariate/Geometric.hs @@ -0,0 +1,82 @@ +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE FunctionalDependencies #-} + +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE KindSignatures #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} + +-- | The method of moments can be used to estimate a number of commonly used distributions. This module is still under construction as I work out the best way to handle morphisms from the Moments3 type to types of other distributions. For more information, see the wikipedia entry: + +module HLearn.Models.Distributions.Univariate.Geometric + ( Geometric + ) + where + +import Control.DeepSeq +import GHC.TypeLits +import qualified Data.Vector.Unboxed as U +import Data.Vector.Unboxed.Deriving + +import qualified Statistics.Distribution as S +import Statistics.Distribution.Geometric + +import HLearn.Algebra +import HLearn.Models.Distributions.Common +import HLearn.Models.Distributions.Univariate.Internal.Moments +import HLearn.Models.Distributions.Visualization.Gnuplot + +------------------------------------------------------------------------------- +-- Geometric + +newtype Geometric sample prob = Geometric { moments :: (Moments3 sample) } + deriving (Read,Show,Eq,Ord,Monoid,Group) + +instance (Num sample) => HomTrainer (Geometric sample prob) where + type Datapoint (Geometric sample prob) = sample + train1dp dp = Geometric $ train1dp dp + +instance (Num sample) => Probabilistic (Geometric sample prob) where + type Probability (Geometric sample prob) = prob + +instance (Integral sample, Floating prob) => PDF (Geometric sample prob) where + pdf dist dp = p*(1-p)^^dp + where + p = geo_p dist + +instance + ( PDF (Geometric sample prob) + , Show prob + , Show sample + , Ord sample + , Ord prob + , Fractional prob + , RealFrac prob + , Integral sample + ) => PlottableDistribution (Geometric sample prob) + where + + plotType _ = Points + + samplePoints dist = [min..max] + where + min = 0 + max = maximum [20,round $ 3*(fromIntegral $ mean dist)] + +geo_p :: (Fractional prob, Integral sample) => Geometric sample prob -> prob +geo_p (Geometric dist) = 1/((fromIntegral $ m1 dist)/(fromIntegral $ m0 dist) +1) + +instance (Integral sample, RealFrac prob) => Mean (Geometric sample prob) where + mean dist = round $ 1/(geo_p dist) + +instance (Integral sample, Fractional prob) => Variance (Geometric sample prob) where + variance dist = (1-p)/(p*p) + where + p = geo_p dist diff --git a/src/HLearn/Models/Distributions/Univariate/Internal/MissingData.hs b/src/HLearn/Models/Distributions/Univariate/Internal/MissingData.hs new file mode 100644 --- /dev/null +++ b/src/HLearn/Models/Distributions/Univariate/Internal/MissingData.hs @@ -0,0 +1,67 @@ +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE FunctionalDependencies #-} +{-# LANGUAGE PolyKinds #-} +{-# LANGUAGE StandaloneDeriving #-} + +-- | Adapts any distribution into one that can handle missing data +module HLearn.Models.Distributions.Univariate.Internal.MissingData + ( MissingData + , MissResponse (..) + ) + where + +import HLearn.Algebra +import HLearn.Models.Distributions.Common + +------------------------------------------------------------------------------- +-- data types + +data MissResponse = Ignore + +newtype MissingData (response :: MissResponse) (basedist:: * -> *) (prob :: *) = + MissingData { basedist :: basedist prob } + deriving (Show,Read,Eq,Ord,Monoid,Group) + +------------------------------------------------------------------------------- +-- Algebra + +instance (Abelian (basedist prob)) => Abelian (MissingData response basedist prob) + +------------------------------------------------------------------------------- +-- Training + +instance + ( HomTrainer (basedist prob) + ) => HomTrainer (MissingData response basedist prob) + where + type Datapoint (MissingData response basedist prob) = Maybe (Datapoint (basedist prob)) + + train1dp Nothing = MissingData mempty + train1dp (Just dp) = MissingData $ train1dp dp + +------------------------------------------------------------------------------- +-- Distribution + +instance Probabilistic (MissingData response basedist prob) where + type Probability (MissingData response basedist prob) = prob + +instance + ( Probability (basedist prob) ~ prob + , PDF (basedist prob) + , Num prob + ) => PDF (MissingData Ignore basedist prob) + where + + {-# INLINE pdf #-} + pdf (MissingData dist) (Just dp) = pdf (dist) dp + pdf (MissingData dist) Nothing = 1 diff --git a/src/HLearn/Models/Distributions/Univariate/Internal/Moments.hs b/src/HLearn/Models/Distributions/Univariate/Internal/Moments.hs new file mode 100644 --- /dev/null +++ b/src/HLearn/Models/Distributions/Univariate/Internal/Moments.hs @@ -0,0 +1,118 @@ +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE FunctionalDependencies #-} + +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE KindSignatures #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} + +-- | The method of moments can be used to estimate a number of commonly used distributions. This module is still under construction as I work out the best way to handle morphisms from the Moments3 type to types of other distributions. For more information, see the wikipedia entry: + +module HLearn.Models.Distributions.Univariate.Internal.Moments + ( Moments3(..) + ) + where + +import Control.DeepSeq +import GHC.TypeLits +import qualified Data.Vector.Unboxed as U +import Data.Vector.Unboxed.Deriving + +import HLearn.Algebra +import HLearn.Models.Distributions.Common +import HLearn.Models.Distributions.Visualization.Gnuplot + +------------------------------------------------------------------------------- +-- Moments3 + +data Moments3 prob = Moments3 + { m0 :: !prob + , m1 :: !prob + , m2 :: !prob + } + deriving (Read,Show,Eq,Ord) + +instance (NFData prob) => NFData (Moments3 prob) where + rnf m = seq m () + +derivingUnbox "Moments3" + [t| (U.Unbox a) => (Moments3 a) -> (a, a, a) |] + [| \ (Moments3 m0 m1 m2) -> (m0,m1,m2) |] + [| \ (m0,m1,m2) -> (Moments3 m0 m1 m2) |] + +------------------------------------------------------------------------------- +-- Algebra + +instance (Num prob) => Abelian (Moments3 prob) +instance (Num prob) => Monoid (Moments3 prob) where + mempty = Moments3 0 0 0 + ma `mappend` mb = Moments3 + { m0 = m0 ma + m0 mb + , m1 = m1 ma + m1 mb + , m2 = m2 ma + m2 mb + } + +instance (Num prob) => Group (Moments3 prob ) where + inverse !m = Moments3 (negate $ m0 m) (negate $ m1 m) (negate $ m2 m) + +instance (Num prob) => HasRing (Moments3 prob) where + type Ring (Moments3 prob) = prob + +instance (Num prob) => Module (Moments3 prob) where + r .* dist = Moments3 + { m0 = r * m0 dist + , m1 = r * m1 dist + , m2 = r * m2 dist + } + + +-- instance (Fractional prob, VU.Unbox prob, SingI n) => LeftModule prob (Moments3 prob n) +-- instance (Fractional prob, VU.Unbox prob) => LeftOperator prob (Moments3 prob n) where +-- (.*) !p !(Moments3 vec) = Moments3 $ VU.map (*p) vec +-- +-- instance (Fractional prob, VU.Unbox prob, SingI n) => RightModule prob (Moments3 prob n) +-- instance (Fractional prob, VU.Unbox prob) => RightOperator prob (Moments3 prob n) where +-- (*.) = flip (.*) +-- +------------------------------------------------------------------------------- +-- Training + +instance (Num prob) => HomTrainer (Moments3 prob) where + type Datapoint (Moments3 prob) = prob + train1dp dp = Moments3 1 dp (dp*dp) + +instance (Num prob) => NumDP (Moments3 prob) where + numdp = m0 + +------------------------------------------------------------------------------- +-- LogNormal + +-- newtype LogNormal prob = LogNormal (Moments3 prob) +-- deriving (Read,Show,Eq,Ord,Semigroup,Monoid,RegularSemigroup) +-- +-- instance ModelParams (LogNormal prob) where +-- type Params (LogNormal prob) = NoParams +-- getparams _ = NoParams +-- +-- instance (Num prob) => HomTrainer (LogNormal prob) where +-- type Datapoint (LogNormal prob) = prob +-- train1dp' params dp = LogNormal $ train1dp' params dp +-- +-- instance (Num prob) => Distribution (LogNormal prob) where +-- type Probability (LogNormal prob) = prob +-- +-- instance (Eq prob, Floating prob) => PDF (LogNormal prob) where +-- pdf (LogNormal dist) 0 = 0 +-- pdf (LogNormal dist) dp = (1/(s*dp*(sqrt $ 2*pi)))*(exp $ (-1)*((log dp)-m)**2/(2*s*s)) +-- where +-- m = 0 +-- s = 1 + + diff --git a/src/HLearn/Models/Distributions/Univariate/LogNormal.hs b/src/HLearn/Models/Distributions/Univariate/LogNormal.hs new file mode 100644 --- /dev/null +++ b/src/HLearn/Models/Distributions/Univariate/LogNormal.hs @@ -0,0 +1,100 @@ +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE FunctionalDependencies #-} + +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE KindSignatures #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} + +-- | LogNormal + +module HLearn.Models.Distributions.Univariate.LogNormal + ( LogNormal + ) + where + +import Debug.Trace + +import HLearn.Algebra +import HLearn.Models.Distributions.Common +import HLearn.Models.Distributions.Univariate.Internal.Moments +import HLearn.Models.Distributions.Visualization.Gnuplot + +------------------------------------------------------------------------------- +-- data types + +newtype LogNormal prob = LogNormal (Moments3 prob) + deriving (Read,Show,Eq,Ord,Monoid,Group) + +------------------------------------------------------------------------------- +-- training + +instance (Floating prob) => HomTrainer (LogNormal prob) where + type Datapoint (LogNormal prob) = prob + train1dp dp = LogNormal $ train1dp $ dp + +------------------------------------------------------------------------------- +-- distribution + +instance Probabilistic (LogNormal prob) where + type Probability (LogNormal prob) = prob + +instance (Floating prob) => PDF (LogNormal prob) where + pdf (LogNormal dist) dp = (1 / (dp * (sqrt $ s2 * 2 * pi)))*(exp $ (-1)*((log dp)-m)^2/(2*s2)) + where +-- sigma2 = variance dist +-- mu = mean dist + + m = 2*log1 - (1/2)*log1 + s2 = log2 - 2*log1 + + log1 = log raw1 + log2 = log raw2 + + raw1 = (m1 dist)/(m0 dist) + raw2 = (m2 dist)/(m0 dist) + +instance (Floating prob) => Mean (LogNormal prob) where + mean (LogNormal dist) = exp $ m+s2/2 + where + m = 2*log1 - (1/2)*log1 + s2 = log2 - 2*log1 + + log1 = log raw1 + log2 = log raw2 + + raw1 = (m1 dist)/(m0 dist) + raw2 = (m2 dist)/(m0 dist) + +instance (Show prob, Floating prob) => Variance (LogNormal prob) where + variance (LogNormal dist) = trace ("m="++show m++"; s2="++show s2) $ ((exp s2) -1)*(exp $ 2*m+s2) + where + m = 2*log1 - (1/2)*log1 + s2 = log2 - 2*log1 + + log1 = log raw1 + log2 = log raw2 + + raw1 = (m1 dist)/(m0 dist) + raw2 = (m2 dist)/(m0 dist) + +instance + ( Floating prob + , Enum prob + , Show prob + , Ord prob + ) => PlottableDistribution (LogNormal prob) where + + plotType _ = Continuous + + samplePoints dist = samplesFromMinMax min max + where + min = (mean dist)-5*(sqrt $ variance dist) + max = (mean dist)+5*(sqrt $ variance dist) diff --git a/src/HLearn/Models/Distributions/Univariate/Normal.hs b/src/HLearn/Models/Distributions/Univariate/Normal.hs new file mode 100644 --- /dev/null +++ b/src/HLearn/Models/Distributions/Univariate/Normal.hs @@ -0,0 +1,81 @@ +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE FunctionalDependencies #-} + +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE KindSignatures #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} + +-- | The method of moments can be used to estimate a number of commonly used distributions. This module is still under construction as I work out the best way to handle morphisms from the Moments3 type to types of other distributions. For more information, see the wikipedia entry: + +module HLearn.Models.Distributions.Univariate.Normal + ( Normal (..) + ) + where + +import Control.DeepSeq +import GHC.TypeLits +import qualified Data.Vector.Unboxed as U +import Data.Vector.Unboxed.Deriving + +import HLearn.Algebra +import HLearn.Models.Distributions.Common +import HLearn.Models.Distributions.Univariate.Internal.Moments +import HLearn.Models.Distributions.Visualization.Gnuplot + +------------------------------------------------------------------------------- +-- data types + +newtype Normal prob = Normal (Moments3 prob) + deriving (Read,Show,Eq,Ord,Monoid,Group,Abelian,Module,NumDP,NFData) + +------------------------------------------------------------------------------- +-- training + +instance (Num prob) => HomTrainer (Normal prob) where + type Datapoint (Normal prob) = prob + train1dp dp = Normal $ train1dp dp + +instance (Num prob) => HasRing (Normal prob) where + type Ring (Normal prob) = prob + +------------------------------------------------------------------------------- +-- algebra + +instance (Num prob) => Probabilistic (Normal prob) where + type Probability (Normal prob) = prob + +instance (Floating prob) => PDF (Normal prob) where + pdf dist dp = (1 / (sqrt $ sigma2 * 2 * pi))*(exp $ (-1)*(dp-mu)*(dp-mu)/(2*sigma2)) + where + sigma2 = variance dist + mu = mean dist + +instance (Fractional prob) => Mean (Normal prob) where + mean (Normal dist) = m1 dist / m0 dist + +instance (Fractional prob) => Variance (Normal prob) where + variance normal@(Normal dist) = m2 dist / m0 dist - (mean normal)*(mean normal) + +instance + ( Floating prob + , Enum prob + , Show prob + , Ord prob + ) => PlottableDistribution (Normal prob) where + + plotType _ = Continuous + + samplePoints dist = samplesFromMinMax min max +-- fmap (\x -> min+x/(numsamples*(max-min))) [0..numsamples] + where + numsamples = 1000 + min = (mean dist)-5*(sqrt $ variance dist) + max = (mean dist)+5*(sqrt $ variance dist) diff --git a/src/HLearn/Models/Distributions/Univariate/Poisson.hs b/src/HLearn/Models/Distributions/Univariate/Poisson.hs new file mode 100644 --- /dev/null +++ b/src/HLearn/Models/Distributions/Univariate/Poisson.hs @@ -0,0 +1,91 @@ +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE FunctionalDependencies #-} + +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE KindSignatures #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} + +-- | The method of moments can be used to estimate a number of commonly used distributions. This module is still under construction as I work out the best way to handle morphisms from the Moments3 type to types of other distributions. For more information, see the wikipedia entry: + +module HLearn.Models.Distributions.Univariate.Poisson + ( Poisson + ) + where + +import Control.DeepSeq +import GHC.TypeLits +import qualified Data.Vector.Unboxed as U +import Data.Vector.Unboxed.Deriving + +import qualified Statistics.Distribution as S +import Statistics.Distribution.Binomial +import Statistics.Distribution.Poisson + +import HLearn.Algebra +import HLearn.Models.Distributions.Common +import HLearn.Models.Distributions.Univariate.Internal.Moments +import HLearn.Models.Distributions.Visualization.Gnuplot + +------------------------------------------------------------------------------- +-- Poisson + +newtype Poisson sample prob = Poisson { pmoments :: (Moments3 sample) } + deriving (Read,Show,Eq,Ord,Monoid,Group) + +instance (Num sample) => HomTrainer (Poisson sample prob) where + type Datapoint (Poisson sample prob) = sample + train1dp dp = Poisson $ train1dp dp + +instance (Num sample) => Probabilistic (Poisson sample prob) where + type Probability (Poisson sample prob) = prob + +-- instance (Integral sample, Floating prob) => PDF (Poisson sample prob) where +-- pdf (Poisson dist) dp +-- | dp < 0 = 0 +-- | dp > 100 = pdf (Normal $ Moments3 (fromIntegral $ m0 dist) (fromIntegral $ m1 dist) (fromIntegral $ m2 dist)) $ fromIntegral dp +-- | otherwise = (lambda^^dp)*(exp $ (-1)*lambda)/(fromIntegral $ factorial dp) +-- where +-- lambda = (fromIntegral $ m1 dist) / (fromIntegral $ m0 dist) +-- factorial :: (Integral a) => a -> a +-- factorial 0 = 1 +-- factorial n = n*(factorial $ n-1) + +instance (Integral sample, Floating prob) => PDF (Poisson sample Double) where + pdf (Poisson dist) dp = S.probability (poisson lambda) $ fromIntegral dp + where + lambda = (fromIntegral $ m1 dist) / (fromIntegral $ m0 dist) + +instance + ( PDF (Poisson sample prob) +-- , PlottableDataPoint sample + , Show prob + , Show sample + , Ord sample + , Ord prob + , Fractional prob + , Integral sample + ) => PlottableDistribution (Poisson sample prob) +-- instance PlottableDistribution (Poisson Int Double) + where + + plotType _ = Points + + samplePoints dist = [min..max] + where + min = 0 + max = maximum [20,floor $ 3*lambda] + lambda = (fromIntegral $ m1 $ pmoments dist) / (fromIntegral $ m0 $ pmoments dist) + +-- instance (Fractional prob) => Mean (Poisson sample prob) where +-- mean (Poisson dist) = (fromIntegral $ m1 $ pmoments dist) / (fromIntegral $ m0 $ pmoments dist) +-- +-- instance (Fractional prob) => Variance (Poisson sample prob) where +-- variance dist = mean dist diff --git a/src/HLearn/Models/Distributions/Visualization/Gnuplot.hs b/src/HLearn/Models/Distributions/Visualization/Gnuplot.hs new file mode 100644 --- /dev/null +++ b/src/HLearn/Models/Distributions/Visualization/Gnuplot.hs @@ -0,0 +1,314 @@ +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE TypeSynonymInstances #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FunctionalDependencies #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE UndecidableInstances #-} + +{-# LANGUAGE OverlappingInstances #-} +-- {-# LANGUAGE IncoherentInstances #-} + +-- | This module contains the functions for plotting distributions using Gnuplot. + +module HLearn.Models.Distributions.Visualization.Gnuplot + ( + -- * Main interface + plotDistribution + + -- ** Plot parameters + , PlotParams(..) + , PicType(..) + , plotFile + , genPlotParams + + -- * Internal + , PlottableDistribution(..) + , Plottable(..) + , PlotType(..) + + -- * Utilities + , samplesFromMinMax + ) + where + +import HLearn.Algebra +import HLearn.Models.Distributions.Common + +import qualified Data.Map as Map +import qualified Data.Vector.Unboxed as VU + +import Control.Monad +import Data.List +import System.IO +import System.Process + +------------------------------------------------------------------------------- +-- data types + +data PlotType = Bar | Points | Continuous + +data PicType = PNG { pngWidth::Int, pngHeight::Int} | EPS + +getExtension :: PicType -> String +getExtension (PNG _ _) = ".png" +getExtension EPS = ".ps" + +data PlotParams = PlotParams + { dataFile :: FilePath + , gnuFile :: FilePath + , picFile :: FilePath + , picType :: PicType + } + +plotFile :: String -> PicType -> PlotParams +plotFile str ext = PlotParams + { dataFile = str++".dat" + , gnuFile = str++".gnu" + , picFile = str++(getExtension ext) + , picType = ext + } + +-- | provided due to backwards compatibility with the nuclear weapons blog post. +genPlotParams :: String -> a -> PlotParams +genPlotParams str a = plotFile str EPS + +------------------------------------------------------------------------------- +-- plotting classes + +class (Show t, Ord t) => Plottable t +instance (Show t, Ord t) => Plottable t + +-- | In order to plot a distribution, it must be an instance of this class. You shouldn't need to know the details. +class + ( Plottable (Datapoint dist) + , Plottable (Probability dist), Num (Probability dist) + , PDF dist + , MaybeShow (Datapoint dist) + ) => PlottableDistribution dist where + + samplePoints :: dist -> [Datapoint dist] + plotType :: dist -> PlotType + + pdfL :: dist -> [Probability dist] + pdfL dist = map (pdf dist) $ samplePoints dist + + plotdata :: dist -> String + plotdata dist = mconcat [maybeshow (x::Datapoint dist) ++ " " ++ show (pdf dist x::Probability dist) ++"\n" | x <- plotPoints] + where + plotPoints = samplePoints dist + +class MaybeShow a where + maybeshow :: a -> String + +instance (Show a) => MaybeShow (Maybe a) where + maybeshow (Just a) = show a + maybeshow Nothing = "Nothing" + +instance (Show a) => MaybeShow a where + maybeshow = show + +------------------------------------------------------------------------------- +-- plotting functions + +-- | Call this function to plot your distribution. You can create the PlotParams manually, or you can use default parameter creating function below. + +plotDistribution :: (PlottableDistribution dist) => PlotParams -> dist -> IO () +plotDistribution params dist = do + -- Create data file + putStrLn "Creating data file..." + datah <- openFile (dataFile params) WriteMode + hPutStrLn datah $ plotdata dist + hClose datah + + -- Create gnuplot file + putStrLn "Plotting data" + gnuh <- openFile (gnuFile params) WriteMode + hPutStrLn gnuh $ gnuplotHeader params + hPutStrLn gnuh $ gnuplot params dist + hClose gnuh + + -- Run gnuplot, creating picture + system $ "gnuplot "++(gnuFile params)++" &" + putStrLn "done." + return () + +gnuplotHeader :: PlotParams -> String +gnuplotHeader params + = terminal + ++ "set output \"" ++ (picFile params) ++ "\" \n" + ++ "unset xtics; unset ytics; unset key \n" +-- ++ "set border 0; set xzeroaxis lt 1; set yzeroaxis lt 1 \n" + ++ "zero(x)=0\n" + ++ "set border 2; set style fill solid 1\n" + where + terminal = case picType params of + EPS -> "set terminal postscript \"Times-Roman\" 25 \n" + ++ "set size 0.81, 1\n" +-- PNG _ _ -> "png" + PNG w h -> "set terminal pngcairo size "++show w++","++show h++" enhanced font 'Times-Roman,10' \n" + +-- class PlotHList t where +-- plotargs :: t -> [(String,String)] -> [String] +-- instance PlotHList (HList '[]) where +-- plotargs HNil _ = [] +-- instance (PlotHList (HList xs)) => PlotHList (HList (x ': xs)) where +-- plotargs (x:::xs) (linecolor,bgcolor):colorL = +-- ["'"++(dataFile params)++"' using 1:"++show i++" lt 1 lw 4 lc rgb '"++linecolor++"' with lines"] +-- :(plotargs xs colorL) + +gnuplot :: (PlottableDistribution dist) => PlotParams -> dist -> String +gnuplot params dist + = yrange + ++ plotcmd + where + plotcmd = case (plotType dist) of + Bar -> "set tics scale 0; set xtics nomirror; unset ytics; unset key\n" + ++ "set xzeroaxis lt 1 lc rgb '#000000'\n" +-- ++ "set ylabel \"Probability\"\n" + ++ "set style data histogram; set style histogram cluster gap 1\n" + ++ "plot '"++(dataFile params)++"' using 2:xticlabels(1) lw 4 linecolor rgb '#0000ff' fs solid 1\n" + Continuous -> "plot '"++(dataFile params)++"' using 1:2:(zero($2)) lt 1 lw 4 lc rgb '#ccccff' with filledcurves, " +-- Continuous -> "plot '"++(dataFile params)++"' using 1:2 lt 1 lw 4 lc rgb '#ccccff' with filledcurves, " + ++ " '"++(dataFile params)++"' using 1:2 lt 1 lw 4 lc rgb '#0000ff' with lines" + Points -> "plot '"++(dataFile params)++"' using 1:2 lt 1 lw 4 lc rgb '#0000ff' with points " + positiveSamples = or $ map (>0) $ pdfL dist + negativeSamples = or $ map (<0) $ pdfL dist + yrange = if positiveSamples && negativeSamples + then "" + else if positiveSamples + then "set yrange [0:]\n" + else "set yrange [:0]\n" + + +-- plotDistributionL :: PlotParams -> [dist] -> IO () +-- plotDistributionL params distL = do +-- -- Create data file +-- putStrLn "Creating data file..." +-- datah <- openFile (dataFile params) WriteMode +-- hPutStrLn datah $ plotdataL distL +-- hClose datah +-- +-- -- Create gnuplot file +-- putStrLn "Plotting data" +-- gnuh <- openFile (gnuFile params) WriteMode +-- hPutStrLn gnuh $ gnuplotL params distL +-- hClose gnuh +-- +-- -- Run gnuplot, creating picture +-- system $ "gnuplot "++(gnuFile params) +-- putStrLn "done." +-- return () +{- + gnuplotL :: PlotParams -> [dist] -> String + gnuplotL params distL + = "set terminal postscript \"Times-Roman\" 25 \n" + ++ "set output \"" ++ (picFile params) ++ "\" \n" + ++ "unset xtics; unset ytics; unset key \n" + ++ "set style fill transparent solid 0.5 noborder\n" + ++ "set border 0; set xzeroaxis lt 1; set yzeroaxis lt 1 \n" + ++ "plot " + -- ++(dataFile params)++"' using 1:2 lt 1 lw 4 lc rgb '#ccccff' with filledcurves, " + ++ (mconcat . intersperse ", " $ + [ -- " '"++(dataFile params)++"' using 1:"++show i++" fs solid 1.0 lc rgb '"++bgcolor++"'" -- "' with filledcurves," + " '"++(dataFile params)++"' using 1:"++show i++" lt 1 lw 4 lc rgb '"++linecolor++"' with lines" + | (i,linecolor,bgcolor,_dist) <- zip4 [2..] linecolorL bgcolorL distL + ]) + where + linecolorL = ["#0000ff", "#00ff00", "#ff0000"] + bgcolorL = ["#ccccff", "#ccffcc", "#ffcccc"]-} +-- plotdataL :: [dist] -> String +-- plotdataL distL = mconcat +-- [show (x::Datapoint dist) ++ (mconcat +-- [ " " ++ show (pdf dist x::Probability dist) +-- | dist <- distL +-- ]) +-- ++"\n" +-- | x <- plotPoints +-- ] +-- where +-- plotPoints = mconcat $ fmap samplePoints distL + + + + +{-instance PlottableDistribution (KDE Double) where + plotdata dist = mconcat [show (x::Double) ++ " " ++ show (pdf dist x::Double) | x <- plotPoints] + where + plotPoints = VU.toList $ samplePoints $ getparams dist + + gnuplot params dist + = "set terminal postscript \"Times-Roman\" 25" + ++ "set output \"" ++ (picFile params) ++ "\"" + ++ "unset xtics; unset ytics; unset key" + ++ "set border 0; set xzeroaxis lt 1; set yzeroaxis lt 1" + ++ "plot '"++(dataFile params)++"' using 1:2 lt 1 lw 4 lc rgb '#ccccff' with filledcurves, \\" + ++ " '"++(dataFile params)++"' using 1:2 lt 1 lw 4 lc rgb '#0000ff' with lines"-} + +{-instance PlottableDistribution (Gaussian Double) where + + minx dist = (meanG dist)-5*(sqrt $ varG dist) + maxx dist = (meanG dist)+5*(sqrt $ varG dist)-} + + +-- instance (Show label, Show prob, Num prob, Ord prob) => PlottableDistribution (Categorical label prob) label where +-- plotdata dist = concat $ do +-- (label,prob) <- dist2list dist +-- return $ show label++" "++show prob++"\n" +-- +-- gnuplot params dist +-- = "set terminal postscript \"Times-Roman\" 25; set output \""++ (picFile params) ++ "\"\n" +-- ++"set tics scale 0; set xtics nomirror; unset ytics; unset key\n" +-- ++"set border 2; set xzeroaxis lt 1\n" +-- ++"set ylabel \"Probability\"\n" +-- ++yrange +-- ++"set style fill border -1\n" +-- ++"set style data histogram; set style histogram cluster gap 1\n" +-- ++"plot '"++(dataFile params)++"' using 2:xticlabels(1)\n" +-- where +-- positiveSamples = or $ map (\(k,v) -> v>0) $ dist2list dist +-- negativeSamples = or $ map (\(k,v) -> v<0) $ dist2list dist +-- yrange = if positiveSamples && negativeSamples +-- then "" +-- else if positiveSamples +-- then "set yrange [0:]\n" +-- else "set yrange [:0]\n" + +{- = "set terminal postscript \"Times-Roman\" 25 \n" + ++ "set output \"" ++ (picFile params) ++ "\" \n" + ++ "unset xtics; unset ytics; unset key\n" + ++ "set border 0; set xzeroaxis lt 1; set yzeroaxis lt 1\n" + ++ "set style data histogram \n" + ++ "plot '"++(dataFile params)++"' using 1:2 \n"-} + +-- 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 () + +samplesFromMinMax min max = fmap (\x -> x/numsamples*(max-min)+min) [0..numsamples] + where + numsamples = 1000 diff --git a/src/HLearn/Models/Distributions/Visualization/Graphviz.hs b/src/HLearn/Models/Distributions/Visualization/Graphviz.hs new file mode 100644 --- /dev/null +++ b/src/HLearn/Models/Distributions/Visualization/Graphviz.hs @@ -0,0 +1,131 @@ +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE FunctionalDependencies #-} +{-# LANGUAGE PolyKinds #-} +{-# LANGUAGE StandaloneDeriving #-} +{-# LANGUAGE OverloadedStrings #-} + +-- | Displays Multivariate dependencies + +module HLearn.Models.Distributions.Visualization.Graphviz + ( MultivariateLabels (..) + , MarkovNetwork (..) + ) + where + +import HLearn.Algebra +import HLearn.Models.Distributions.Multivariate.Interface +import HLearn.Models.Distributions.Multivariate.Internal.CatContainer +import HLearn.Models.Distributions.Multivariate.Internal.Container +import HLearn.Models.Distributions.Multivariate.Internal.TypeLens + +import Data.GraphViz.Exception +import Data.GraphViz hiding (graphToDot) +import Data.GraphViz.Attributes.Complete{-( Attribute(RankDir, Splines, FontName) + , RankDir(FromLeft), EdgeType(SplineEdges))-} +import Control.Arrow(second) +import GHC.TypeLits + +------------------------------------------------------------------------------- +-- clases + +class (Trainable datatype) => MultivariateLabels datatype where + getLabels :: datatype -> [String] + +class (MultivariateLabels (Datapoint dist)) => MarkovNetwork dist where + graphL :: dist -> [String] -> [(String,[String])] + + plotNetwork :: FilePath -> dist -> IO Bool + plotNetwork file dist = graphToDotPng file $ graphL dist $ getLabels (undefined :: Datapoint dist) + +------------------------------------------------------------------------------- +-- instances + +instance + ( MultivariateLabels datapoint + ) => MarkovNetwork (Multivariate datapoint '[] prob) + where + graphL _ labels = [] + +instance + ( MultivariateLabels datapoint + , MarkovNetwork (Multivariate datapoint xs prob) + ) => MarkovNetwork (Multivariate datapoint ( ('[]) ': xs) prob) + where + graphL _ labels = graphL (undefined :: Multivariate datapoint xs prob) labels + +instance + ( MultivariateLabels datapoint + , MarkovNetwork (Multivariate datapoint ( ys ': xs) prob) + ) => MarkovNetwork (Multivariate datapoint ( (Ignore' label ': ys) ': xs) prob) + where + graphL _ labels = (graphL (undefined :: Multivariate datapoint ( ys ': xs) prob) (tail labels)) + +instance + ( MultivariateLabels datapoint + , MarkovNetwork (Multivariate datapoint ( ys ': xs) prob) + ) => MarkovNetwork (Multivariate datapoint ( (CatContainer label ': ys) ': xs) prob) + where + graphL _ labels = (head labels, tail labels) + : (graphL (undefined :: Multivariate datapoint ( ys ': xs) prob) (tail labels)) + +instance + ( MultivariateLabels datapoint + , MarkovNetwork (Multivariate datapoint (ys ': xs) prob) + ) => MarkovNetwork (Multivariate datapoint ( (Container dist label ': ys) ': xs) prob) + where + graphL _ l = (head l,[]):(graphL (undefined::Multivariate datapoint (ys ': xs) prob) (tail l)) + +instance + ( MultivariateLabels datapoint + , SingI (Length labelL) + , MarkovNetwork (Multivariate datapoint ( ys ': xs) prob) + ) => MarkovNetwork (Multivariate datapoint ( (MultiContainer dist (labelL:: [*]) ': ys) ': xs) prob) + where + graphL _ l = go (take n l) ++ (graphL (undefined :: Multivariate datapoint ( ys ': xs ) prob) $ drop n l) + where + go [] = [] + go (x:xs) = (x,xs):(go xs) + + n = fromIntegral $ fromSing $ (sing :: Sing (Length labelL)) + +------------------------------------------------------------------------------- +-- Graphviz helpers + +--------------------------------------- +-- These functions are taken from the graphviz tutorial at: +-- http://ivanmiljenovic.wordpress.com/2011/10/16/graphviz-in-vacuum/ + +graphToDot :: (Ord a) => [(a, [a])] -> DotGraph a +graphToDot = graphToDotParams vacuumParams + +graphToDotParams :: (Ord a, Ord cl) => GraphvizParams a () () cl l -> [(a, [a])] -> DotGraph a +graphToDotParams params nes = graphElemsToDot params ns es + where + ns = map (second $ const ()) nes + es = concatMap mkEs nes + mkEs (f,ts) = map (\t -> (f,t,())) ts + +------------------------------------------------ + +vacuumParams :: GraphvizParams a () () () () +vacuumParams = defaultParams { globalAttributes = gStyle } + +gStyle :: [GlobalAttributes] +gStyle = [ GraphAttrs [RankDir FromLeft, {-Splines SplineEdges, -}FontName "courier", Layout Circo] + , NodeAttrs [textLabel "\\N", shape PlainText, fontColor Black, Shape Ellipse, style filled, fillColor AliceBlue, penWidth 2, color Navy] + , EdgeAttrs [color Black, Dir NoDir] + ] + +graphToDotPng :: FilePath -> [(String,[String])] -> IO Bool +graphToDotPng fpre g = handle (\(e::GraphvizException) -> return False) + $ addExtension (runGraphviz (graphToDot g)) Png fpre >> return True diff --git a/src/examples/Criterion.hs b/src/examples/Criterion.hs deleted file mode 100644 --- a/src/examples/Criterion.hs +++ /dev/null @@ -1,66 +0,0 @@ -import Criterion.Main - -import Statistics.Distribution.Normal -import qualified Data.Vector as V -import qualified Data.Vector.Unboxed as VU -import qualified Data.Vector.Generic as G - -import HLearn.Algebra -import HLearn.Models.Distributions --- import qualified HLearn.Models.Distributions.GaussianOld as GO --- import qualified HLearn.Models.Distributions.GaussianOld2 as GO2 - -import qualified Control.ConstraintKinds as CK - -size = 10^8 -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)) --- , bench "batch train VU 1e6" $ nf ((batch train) GaussianParams) (VU.enumFromN (0::Double) [0..1e6::Double]) --- , bench "parallel2 batch train [] 1e6" $ nf ((parallel $ batch train) GaussianParams) (V.enumFromN (0::Double) (10^7)) - ,-} --- bench "Parallel GaussianOld2" $ nf (parallel $ batch (train GO2.GaussianParams)) (VU.enumFromN (0::Double) (10^8)) --- , bench "GaussianOld2" $ nf (batch (train GO2.GaussianParams)) (VU.enumFromN (0::Double) (10^8)) --- bench "Parallel GaussianOld" $ nf (parallel $ batch (train GO.GaussianParams)) (VU.enumFromN (0::Double) (10^8)) --- , bench "GaussianOld" $ nf (batch (train GO.GaussianParams)) (VU.enumFromN (0::Double) (10^8)) --- bench "Parallel Gaussian2" $ nf (parallel $ batch (train GaussianParams)) (VU.enumFromN (0::Double) (10^8)) --- , bench "Gaussian2" $ nf (batch (train GaussianParams)) (VU.enumFromN (0::Double) (10^8)) - - {-bench "Parallel Gaussian" $ nf (parallel $ batch (trainSG)) (VU.enumFromN (0::Double) (10^8)) - , bench "Gaussian" $ nf (batch (trainSG)) (VU.enumFromN (0::Double) (10^8)) - ,-} -- bench "Categorical" $ nf (train :: [String] -> Categorical String Double) (concat $ replicate (10^4) ["a","b"]) --- , bench "normalFromSample - 1e5" $ whnf (normalFromSample . VU.enumFromN 1) (10^8) - - --- instance NFData NormalDistribution where --- rnf (NormalDistribution a b c d) = deepseq a $ deepseq b $ deepseq c $ rnf d - --- parallel3 strat train = \modelparams datapoint -> --- foldl1 (<>) $ map (CK.foldl1 (<>) . CK.fmap (train modelparams . CK.pure)) (CK.partition 2 datapoint) - - --- test (ds1,ds2) = m1 <> m2 --- where --- [m1,m2] = parMap rdeepseq ((batch train) GaussianParams) [ds1,ds2] --- -- [m1,m2] = parMap rdeepseq ((batch train) GaussianParams) [ds1,ds2] --- --- -- test (ds1,ds2) = runEval $ do --- -- let [m1,m2] = {-parMap rdeepseq-} map (batch train GaussianParams) [ds1,ds2] --- -- -- m1 <- rparWith rseq $ batch train GaussianParams ds1 --- -- -- m2 <- rparWith rseq $ batch train GaussianParams ds2 --- -- {- let m1 = batch train GaussianParams ds1 --- -- let m2 = batch train GaussianParams ds2-} --- -- return $ m1 <> m2 --- --- dsL = replicate 4 ds --- -- dsL = [ds1,ds2] --- ds = [0..1000000::Double] --- ds1 = V.enumFromN (0::Double) (10^7) --- ds2 = V.enumFromN (0.5::Double) (10^7) diff --git a/src/examples/SpaceTests.hs b/src/examples/SpaceTests.hs deleted file mode 100644 --- a/src/examples/SpaceTests.hs +++ /dev/null @@ -1,24 +0,0 @@ -import HLearn.Algebra -import HLearn.Models.Distributions - -import Control.Monad -import Control.Monad.ST -import Data.List -import Data.STRef -import Data.Number.LogFloat - -import Statistics.Distribution -import Statistics.Distribution.Normal - -import qualified Data.Vector as V -import qualified Data.Vector.Unboxed as VU -import qualified Data.Vector.Generic as G -import qualified Data.Foldable as F - -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 $ normalFromSample $ VU.enumFromN 0 size - putStrLn "done."