diff --git a/Goal/Probability.hs b/Goal/Probability.hs --- a/Goal/Probability.hs +++ b/Goal/Probability.hs @@ -1,13 +1,40 @@ +{-# OPTIONS_GHC -fplugin=GHC.TypeLits.KnownNat.Solver -fplugin=GHC.TypeLits.Normalise -fconstraint-solver-iterations=10 #-} +{-# LANGUAGE + RankNTypes, + TypeOperators, + FlexibleContexts, + ScopedTypeVariables +#-} +-- | The main module of goal-probability. Import this module to use all the +-- types, functions, and classes provided by goal-probability. module Goal.Probability - ( module System.Random.MWC - , module System.Random.MWC.Monad - , module Goal.Probability.Statistical + ( -- * Package Exports + module Goal.Probability.Statistical , module Goal.Probability.ExponentialFamily + , module Goal.Probability.Conditional , module Goal.Probability.Distributions - , module Goal.Probability.Graphical - , module Goal.Probability.Graphical.Harmonium - , module Goal.Probability.Graphical.NeuralNetwork - , module Goal.Probability + , module Goal.Probability.Distributions.Gaussian + , module Goal.Probability.Distributions.CoMPoisson + -- * Stochastic Operations + , shuffleList + , resampleVector + , subsampleVector + , noisyFunction + -- ** Circuits + , minibatcher + -- * Statistics + , estimateMeanVariance + , estimateMeanSD + , estimateFanoFactor + , estimateCoefficientOfVariation + , estimateCorrelation + , estimateCorrelations + , histograms + -- ** Model Selection + , akaikesInformationCriterion + , bayesianInformationCriterion + --, conditionalAkaikesInformationCriterion + --, conditionalBayesianInformationCriterion ) where @@ -16,58 +43,191 @@ -- Re-exports -- -import System.Random.MWC hiding (uniform,uniformR) -import System.Random.MWC.Monad hiding (save) - -import qualified System.Random.MWC.Monad as S (save) - import Goal.Probability.Statistical import Goal.Probability.ExponentialFamily +import Goal.Probability.Conditional import Goal.Probability.Distributions -import Goal.Probability.Graphical -import Goal.Probability.Graphical.Harmonium -import Goal.Probability.Graphical.NeuralNetwork +import Goal.Probability.Distributions.Gaussian +import Goal.Probability.Distributions.CoMPoisson -- Package -- import Goal.Core import Goal.Geometry +import qualified Goal.Core.Vector.Boxed as B +import qualified Goal.Core.Vector.Storable as S +import qualified Goal.Core.Vector.Generic.Mutable as M +import qualified Goal.Core.Vector.Generic as G +import qualified Data.Vector.Generic.Mutable.Base as MV +import qualified Data.Vector as V + +import qualified Statistics.Sample as STAT hiding (range) +import qualified Statistics.Sample.Histogram as STAT +import qualified Data.Vector.Storable as VS + +import qualified System.Random.MWC as R +import qualified System.Random.MWC.Distributions as R + + +--- Statistics --- + + +-- | Estimate the mean and variance of a sample (with Bessel's correction) +estimateMeanVariance + :: Traversable f + => f Double + -> (Double,Double) +estimateMeanVariance xs = STAT.meanVarianceUnb . VS.fromList $ toList xs + +-- | Estimate the mean and variance of a sample (with Bessel's correction) +estimateMeanSD + :: Traversable f + => f Double + -> (Double,Double) +estimateMeanSD xs = + let (mu,vr) = estimateMeanVariance xs + in (mu,sqrt vr) + +-- | Estimate the Fano Factor of a sample. +estimateFanoFactor + :: Traversable f + => f Double + -> Double +estimateFanoFactor xs = + let (mu,vr) = estimateMeanVariance xs + in vr / mu + +-- | Estimate the coefficient of variation from a sample. +estimateCoefficientOfVariation :: Traversable f => f Double -> Double +estimateCoefficientOfVariation zs = + let (mu,vr) = estimateMeanVariance zs + in sqrt vr / mu + +-- | Computes the empirical covariance matrix given a sample if iid random vectors. +estimateCorrelations + :: forall k x v . (G.VectorClass v x, G.VectorClass v Double, KnownNat k, Real x) + => [G.Vector v k x] + -> S.Matrix k k Double +estimateCorrelations zs = + let mnrm :: Source # MultivariateNormal k + mnrm = mle $ G.convert . G.map realToFrac <$> zs + in multivariateNormalCorrelations mnrm + +-- | Computes the empirical covariance matrix given a sample from a bivariate random variable. +estimateCorrelation + :: [(Double,Double)] + -> Double +estimateCorrelation zs = STAT.correlation $ V.fromList zs + +-- | Computes histograms (and densities) with the given number of bins for the +-- given list of samples. Bounds can be given or computed automatically. The +-- returned values are the list of bin centres and the binned samples. If bounds +-- are given but are not greater than all given sample points, then an error +-- will be thrown. +histograms + :: Int -- ^ Number of Bins + -> Maybe (Double, Double) -- ^ Maybe bin bounds + -> [[Double]] -- ^ Datasets + -> ([Double],[[Int]],[[Double]]) -- ^ Bin centres, counts, and densities for each dataset +histograms nbns mmnmx smpss = + let (mn,mx) = case mmnmx of + Just (mn0,mx0) -> (mn0,mx0) + Nothing -> STAT.range nbns . VS.fromList $ concat smpss + stp = (mx - mn) / fromIntegral nbns + bns = take nbns [ mn + stp/2 + stp * fromIntegral n | n <- [0 :: Int,1..] ] + hsts = VS.toList . STAT.histogram_ nbns mn mx . VS.fromList <$> smpss + ttls = sum <$> hsts + dnss = do + (hst,ttl) <- zip hsts ttls + return $ if ttl == 0 + then [] + else (/(fromIntegral ttl * stp)) . fromIntegral <$> hst + in (bns,hsts,dnss) + + --- Stochastic Functions --- -seed :: RandST s Seed --- | This little guy creates a seed. It's necessary to avoid name space --- collisions. -seed = S.save +-- | Shuffle the elements of a list. +shuffleList :: [a] -> Random [a] +shuffleList xs = V.toList <$> Random (R.uniformShuffle (V.fromList xs)) -randomElement :: [x] -> RandST r x --- | Returns a random element from a list. -randomElement xs = do - u <- uniform - let elm = round $ fromIntegral (length xs - 1) * (u :: Double) - return $ xs !! elm +-- | A 'Circuit' that helps fitting data based on minibatches. Essentially, it +-- creates an infinite list out of shuffled versions of the input list, and +-- breaks down and returns the result in chunks of the specified size. +minibatcher :: Int -> [x] -> Chain Random [x] +minibatcher nbtch xs0 = accumulateFunction [] $ \() xs -> + if length (take nbtch xs) < nbtch + then do + xs1 <- shuffleList xs0 + let (hds',tls') = splitAt nbtch (xs ++ xs1) + return (hds',tls') + else do + let (hds',tls') = splitAt nbtch xs + return (hds',tls') -noisyFunction :: (Generative c m, Num (Sample m)) - => (c :#: m) -- ^ Noise model - -> (x -> Sample m) -- ^ Function - -> x - -> RandST r (Sample m) +-- | Returns a uniform sample of elements from the given vector with replacement. +resampleVector :: (KnownNat n, KnownNat k) => B.Vector n x -> Random (B.Vector k x) +resampleVector xs = do + ks <- B.replicateM $ Random (R.uniformR (0, B.length xs-1)) + return $ B.backpermute xs ks + -- | Returns a sample from the given function with added noise. +noisyFunction + :: (Generative c x, Num (SamplePoint x)) + => Point c x -- ^ Noise model + -> (y -> SamplePoint x) -- ^ Function + -> y -- ^ Input + -> Random (SamplePoint x) -- ^ Stochastic Output noisyFunction m f x = do - ns <- generate m + ns <- samplePoint m return $ f x + ns -noisyRange - :: Double -- ^ The min of the function input - -> Double -- ^ The max function input - -> Int -- ^ Number of samples to draw from the function - -> Double -- ^ Standard deviation of the noise - -> (Double -> Double) -- ^ Mixture function - -> RandST s [(Double,Double)] -{-| Returns a set of samples from the given function with additive Gaussian noise. -} -noisyRange mn mx n sd f = do - let xs = range mn mx n - d = chart Standard $ fromList Normal [0,sd^2] - fxs <- mapM (\x -> (+ f x) <$> generate d) xs - return $ zip xs fxs +-- | Take a random, unordered subset of a list. +subsampleVector + :: forall k m v x . (KnownNat k, KnownNat m, G.VectorClass v x) + => G.Vector v (k + m) x + -> Random (G.Vector v k x) +subsampleVector v = Random $ \gn -> do + let k = natValInt (Proxy :: Proxy k) + mv <- G.thaw v + randomSubSample0 k mv gn + v' <- G.unsafeFreeze mv + let foo :: (G.Vector v k x, G.Vector v m x) + foo = G.splitAt v' + return $ fst foo + +randomSubSample0 + :: (KnownNat n, PrimMonad m, MV.MVector v a) + => Int -> G.MVector v n (PrimState m) a -> R.Gen (PrimState m) -> m () +randomSubSample0 k v gn = looper 0 + where n = M.length v + looper i + | i == k = return () + | otherwise = do + j <- R.uniformR (i,n-1) gn + M.unsafeSwap v i j + looper (i+1) + + +-- | Calculate the AIC for a given model and sample. +akaikesInformationCriterion + :: forall c x s . (Manifold x, LogLikelihood c x s) + => c # x + -> [s] + -> Double +akaikesInformationCriterion p xs = + let d = natVal (Proxy :: Proxy (Dimension x)) + in 2 * fromIntegral d - 2 * logLikelihood xs p * fromIntegral (length xs) + +-- | Calculate the BIC for a given model and sample. +bayesianInformationCriterion + :: forall c x s . (LogLikelihood c x s, Manifold x) + => c # x + -> [s] + -> Double +bayesianInformationCriterion p xs = + let d = natVal (Proxy :: Proxy (Dimension x)) + n = length xs + in log (fromIntegral n) * fromIntegral d - 2 * logLikelihood xs p * fromIntegral (length xs) diff --git a/Goal/Probability/Conditional.hs b/Goal/Probability/Conditional.hs new file mode 100644 --- /dev/null +++ b/Goal/Probability/Conditional.hs @@ -0,0 +1,228 @@ +{-# OPTIONS_GHC -fplugin=GHC.TypeLits.KnownNat.Solver -fplugin=GHC.TypeLits.Normalise -fconstraint-solver-iterations=10 #-} +{-# LANGUAGE UndecidableInstances #-} + +-- | 'Statistical' models where the observations depend on known conditions. +module Goal.Probability.Conditional + ( SampleMap + -- ** Markov Kernels + , (>.>*) + , (>$>*) + , (*<.<) + , (*<$<) + -- ** Conditional Distributions + , conditionalLogLikelihood + , conditionalLogLikelihoodDifferential + , conditionalDataMap + , kFoldMap + , kFoldMap' + --, mapToConditionalData + , mapConditionalLogLikelihood + , mapConditionalLogLikelihoodDifferential + , parMapConditionalLogLikelihood + , parMapConditionalLogLikelihoodDifferential + ) where + + +--- Imports --- + + +-- Goal -- + +import Goal.Core +import Goal.Geometry + +import Goal.Probability.Statistical +import Goal.Probability.ExponentialFamily + +import qualified Data.Map.Strict as M +import qualified Data.List as L + +import Control.Parallel.Strategies + + +--- Generic --- + + +-- | Evalutes the given conditional distribution at a 'SamplePoint'. +(>.>*) :: (Map Natural f y x, ExponentialFamily x) + => Natural # f y x + -> SamplePoint x + -> Natural # y +(>.>*) p x = p >.> sufficientStatistic x + +-- | Mapped application of conditional distributions on a 'Sample'. +(>$>*) :: (Map Natural f y x, ExponentialFamily x) + => Natural # f y x + -> Sample x + -> [Natural # y] +(>$>*) p xs = p >$> (sufficientStatistic <$> xs) + +infix 8 >.>* +infix 8 >$>* + +-- | Applies the transpose of a 'Bilinear' 'Map' to a 'SamplePoint'. +(*<.<) :: (Map Natural f x y, Bilinear f y x, ExponentialFamily y) + => SamplePoint y + -> Natural # f y x + -> Natural # x +(*<.<) x p = sufficientStatistic x <.< p + +-- | Mapped transpose application on a 'Sample'. +(*<$<) :: (Map Natural f x y, Bilinear f y x, ExponentialFamily y) + => Sample y + -> Natural # f y x + -> [Natural # x] +(*<$<) xs p = (sufficientStatistic <$> xs) <$< p + +infix 8 *<.< +infix 8 *<$< + + +-- | A synonym for Maps from Inputs to Outputs that matches the confusing, +-- backwards style of Goal. +type SampleMap z x = M.Map (SamplePoint x) (Sample z) + + +dependantLogLikelihood + :: (LogLikelihood Natural y s, Map Natural f y x) + => [([s], Mean # x)] -> Natural # f y x -> Double +dependantLogLikelihood ysxs chrm = + let (yss,xs) = unzip ysxs + in average . zipWith logLikelihood yss $ chrm >$> xs + +dependantLogLikelihoodDifferential + :: (LogLikelihood Natural y s, Propagate Natural f y x) + => [([s], Mean # x)] -> Natural # f y x -> Mean # f y x +dependantLogLikelihoodDifferential ysxs chrm = + let (yss,xs) = unzip ysxs + (df,yhts) = propagate mys xs chrm + mys = zipWith logLikelihoodDifferential yss yhts + in df + +dependantLogLikelihoodPar + :: (LogLikelihood Natural y s, Map Natural f y x) + => [([s], Mean # x)] -> Natural # f y x -> Double +dependantLogLikelihoodPar ysxs chrm = + let (yss,xs) = unzip ysxs + in average . parMap rdeepseq (uncurry logLikelihood) . zip yss $ chrm >$> xs + +dependantLogLikelihoodDifferentialPar + :: (LogLikelihood Natural y s, Propagate Natural f y x) + => [([s], Mean # x)] -> Natural # f y x -> Mean # f y x +dependantLogLikelihoodDifferentialPar ysxs chrm = + let (yss,xs) = unzip ysxs + (df,yhts) = propagate mys xs chrm + mys = parMap rdeepseq (uncurry logLikelihoodDifferential) $ zip yss yhts + in df + +-- | Turns a list of input/output pairs into a Map, by collecting into lists the +-- different outputs to each particular input. +conditionalDataMap + :: Ord x + => [(t, x)] -- ^ Output/Input Pairs + -> M.Map x [t] -- ^ Input Output map +conditionalDataMap = foldl' folder M.empty + where folder mp (t,x) = + let ts = M.lookup x mp + ts' = maybe [t] (t:) ts + in M.insert x ts' mp + --M.fromListWith (++) [(x, [y]) | (y, x) <- yxs] + +-- | Partition a conditional dataset into k > 1 (training,validation) pairs, +-- where each dataset condition is partitioned to match its size. +kFoldMap + :: Ord x => Int -> M.Map x [y] -> [(M.Map x [y], M.Map x [y])] +kFoldMap k ixzmp = + let ixzmps = kFold k <$> ixzmp + ixs = M.keys ixzmp + tvzss = M.elems ixzmps + tvxzmps = M.fromList . zip ixs <$> L.transpose tvzss + in zip (fmap fst <$> tvxzmps) (fmap snd <$> tvxzmps) + +-- | Partition a conditional dataset into k > 2 (training,test,validation) triplets, +-- where each dataset condition is partitioned to match its size. +kFoldMap' + :: Ord x => Int -> M.Map x [y] -> [(M.Map x [y], M.Map x [y], M.Map x [y])] +kFoldMap' k ixzmp = + let ixzmps = kFold' k <$> ixzmp + ixs = M.keys ixzmp + tvzss = M.elems ixzmps + tvxzmps = M.fromList . zip ixs <$> L.transpose tvzss + in zip3 (fmap (\(x,_,_) -> x) <$> tvxzmps) + (fmap (\(_,x,_) -> x) <$> tvxzmps) + (fmap (\(_,_,x) -> x) <$> tvxzmps) + +--mapToConditionalData :: M.Map x [y] -> [(y,x)] +--mapToConditionalData mp = +-- let (xs,zss) = unzip $ M.toAscList mp +-- in concat $ zipWith (\x zs -> zip zs $ repeat x) xs zss + + +-- | The conditional 'logLikelihood' for a conditional distribution. +conditionalLogLikelihood + :: (ExponentialFamily x, Map Natural f y x, LogLikelihood Natural y t) + => [(t, SamplePoint x)] -- ^ Output/Input Pairs + -> Natural # f y x -- ^ Function + -> Double -- ^ conditional cross entropy estimate +conditionalLogLikelihood yxs f = + let ysxs = [ ([y],sufficientStatistic x) | (y,x) <- yxs ] + in dependantLogLikelihood ysxs f + +-- | The conditional 'logLikelihoodDifferential' for a conditional distribution. +conditionalLogLikelihoodDifferential + :: ( ExponentialFamily x, LogLikelihood Natural y t, Propagate Natural f y x ) + => [(t, SamplePoint x)] -- ^ Output/Input Pairs + -> Natural # f y x -- ^ Function + -> Mean # f y x -- ^ Differential +conditionalLogLikelihoodDifferential yxs f = + let ysxs = [ ([y],sufficientStatistic x) | (y,x) <- yxs ] + in dependantLogLikelihoodDifferential ysxs f + +-- | The conditional 'logLikelihood' for a conditional distribution, where +-- redundant conditions/inputs are combined. This can dramatically increase performance when +-- the number of distinct conditions/inputs is small. +mapConditionalLogLikelihood + :: ( ExponentialFamily x, Map Natural f y x, LogLikelihood Natural y t ) + => M.Map (SamplePoint x) [t] -- ^ Output/Input Pairs + -> Natural # f y x -- ^ Function + -> Double -- ^ conditional cross entropy estimate +mapConditionalLogLikelihood xtsmp = + dependantLogLikelihood [ (ts, sufficientStatistic x) | (x,ts) <- M.toList xtsmp] + +-- | The conditional 'logLikelihoodDifferential', where redundant conditions are +-- combined. This can dramatically increase performance when the number of +-- distinct conditions is small. +mapConditionalLogLikelihoodDifferential + :: ( ExponentialFamily x, LogLikelihood Natural y t + , Propagate Natural f y x, Ord (SamplePoint x) ) + => M.Map (SamplePoint x) [t] -- ^ Output/Input Pairs + -> Natural # f y x -- ^ Function + -> Mean # f y x -- ^ Differential +mapConditionalLogLikelihoodDifferential xtsmp = + dependantLogLikelihoodDifferential [ (ts, sufficientStatistic x) | (x,ts) <- M.toList xtsmp] + +-- | The conditional 'logLikelihood' for a conditional distribution, where +-- redundant conditions/inputs are combined. This can dramatically increase performance when +-- the number of distinct conditions/inputs is small. +parMapConditionalLogLikelihood + :: ( ExponentialFamily x, Map Natural f y x, LogLikelihood Natural y t ) + => M.Map (SamplePoint x) [t] -- ^ Output/Input Pairs + -> Natural # f y x -- ^ Function + -> Double -- ^ conditional cross entropy estimate +parMapConditionalLogLikelihood xtsmp = + dependantLogLikelihoodPar [ (ts, sufficientStatistic x) | (x,ts) <- M.toList xtsmp] + +-- | The conditional 'logLikelihoodDifferential', where redundant conditions are +-- combined. This can dramatically increase performance when the number of +-- distinct conditions is small. +parMapConditionalLogLikelihoodDifferential + :: ( ExponentialFamily x, LogLikelihood Natural y t + , Propagate Natural f y x, Ord (SamplePoint x) ) + => M.Map (SamplePoint x) [t] -- ^ Output/Input Pairs + -> Natural # f y x -- ^ Function + -> Mean # f y x -- ^ Differential +parMapConditionalLogLikelihoodDifferential xtsmp = + dependantLogLikelihoodDifferentialPar [ (ts, sufficientStatistic x) | (x,ts) <- M.toList xtsmp] + + + diff --git a/Goal/Probability/Distributions.hs b/Goal/Probability/Distributions.hs --- a/Goal/Probability/Distributions.hs +++ b/Goal/Probability/Distributions.hs @@ -1,18 +1,21 @@ --- | Various instances of 'Statistical' 'Manifold's. -module Goal.Probability.Distributions ( - -- * General Statistical Manifolds - CurvedCategorical (CurvedCategorical) - , Uniform (Uniform) - -- * Exponential Family Manifolds - , Bernoulli (Bernoulli) - , Binomial (Binomial) - , Categorical (Categorical) - , Poisson (Poisson) - , Normal (Normal) - , MeanNormal (MeanNormal) - , MultivariateNormal (MultivariateNormal) - -- * Util - , muSigmaToMultivariateNormal +{-# OPTIONS_GHC -fplugin=GHC.TypeLits.KnownNat.Solver -fplugin=GHC.TypeLits.Normalise -fconstraint-solver-iterations=10 #-} +{-# LANGUAGE UndecidableInstances,TypeApplications #-} + +-- | Various instances of statistical manifolds, with a focus on exponential +-- families. In the documentation we use \(X\) to indicate a random variable +-- with the distribution being documented. +module Goal.Probability.Distributions + ( -- * Univariate + Bernoulli + , Binomial + , Categorical + , categoricalWeights + , Poisson + , VonMises + -- * Multivariate + , Dirichlet + -- * LocationShape + , LocationShape (LocationShape) ) where -- Package -- @@ -23,628 +26,588 @@ import Goal.Geometry --- Qualified -- +import qualified Goal.Core.Vector.Storable as S +import qualified Goal.Core.Vector.Boxed as B +import qualified Goal.Core.Vector.Generic as G -import qualified Data.Vector.Storable as C -import qualified Numeric.LinearAlgebra.HMatrix as M +import qualified Numeric.GSL.Special.Bessel as GSL +import qualified Numeric.GSL.Special.Gamma as GSL +import qualified Numeric.GSL.Special.Psi as GSL +import qualified System.Random.MWC as R +import qualified System.Random.MWC.Distributions as R --- Unqualified -- +import Foreign.Storable -import System.Random.MWC.Monad -import System.Random.MWC.Distributions.Monad -import Statistics.Sample hiding (mean) -import Numeric.SpecFunctions +-- Location Shape -- +-- | A 'LocationShape' 'Manifold' is a 'Product' of some location 'Manifold' and +-- some shape 'Manifold'. +newtype LocationShape l s = LocationShape (l,s) + +deriving instance (Manifold l, Manifold s) => Manifold (LocationShape l s) +deriving instance (Manifold l, Manifold s) => Product (LocationShape l s) + -- Uniform -- -data Uniform = Uniform Double Double deriving (Eq, Read, Show) +-- Bernoulli Distribution -- -instance Manifold Uniform where - dimension _ = 0 +-- | The Bernoulli family with 'Bool'ean 'SamplePoint's. (because why not). The source coordinate is \(P(X = True)\). +data Bernoulli -instance Statistical Uniform where - type SampleSpace Uniform = Continuum - sampleSpace _ = Continuum +-- Binomial Distribution -- -instance Generative Standard Uniform where - generate p = - let (Uniform a b) = manifold p - in uniformR (a,b) +-- | A distribution over the sum of 'True' realizations of @n@ 'Bernoulli' +-- random variables. The 'Source' coordinate is the probability of \(P(X = True)\) +-- for each 'Bernoulli' random variable. +data Binomial (n :: Nat) -instance AbsolutelyContinuous Standard Uniform where - density p x = - let (Uniform a b) = manifold p - in if x >= a && x <= b - then recip $ b - a - else 0 +-- | Returns the number of trials used to define this binomial distribution. +binomialTrials :: forall c n. KnownNat n => Point c (Binomial n) -> Int +binomialTrials _ = natValInt (Proxy :: Proxy n) --- Bernoulli Distribution -- +-- | Returns the number of trials used to define this binomial distribution. +binomialSampleSpace :: forall n . KnownNat n => Proxy (Binomial n) -> Int +binomialSampleSpace _ = natValInt (Proxy :: Proxy n) --- | The Bernoulli 'Family' with 'SampleSpace' 'Bernoulli' = 'Bool' (because why not). -data Bernoulli = Bernoulli deriving (Eq, Read, Show) +-- Categorical Distribution -- -instance Manifold Bernoulli where - dimension _ = 1 +-- | A 'Categorical' distribution where the probability of the first category +-- \(P(X = 0)\) is given by the normalization constraint. +data Categorical (n :: Nat) -instance Statistical Bernoulli where - type SampleSpace Bernoulli = Boolean - sampleSpace Bernoulli = Boolean +-- | Takes a weighted list of elements representing a probability mass function, and +-- returns a sample from the Categorical distribution. +sampleCategorical :: KnownNat n => S.Vector n Double -> Random Int +sampleCategorical ps = do + let ps' = S.postscanl' (+) 0 ps + p <- Random R.uniform + let midx = (+1) . finiteInt <$> S.findIndex (> p) ps' + return $ fromMaybe 0 midx -instance Generative Standard Bernoulli where - generate p = bernoulli . C.head $ coordinates p +-- | Returns the probabilities over the whole sample space \((0 \ldots n)\) of the +-- given categorical distribution. +categoricalWeights + :: Transition c Source (Categorical n) + => c # Categorical n + -> S.Vector (n+1) Double +categoricalWeights wghts0 = + let wghts = coordinates $ toSource wghts0 + in S.cons (1-S.sum wghts) wghts -instance AbsolutelyContinuous Standard Bernoulli where - density p True = C.head $ coordinates p - density p False = 1 - C.head (coordinates p) +-- | A 'Dirichlet' manifold contains distributions over weights of a +-- 'Categorical' distribution. +data Dirichlet (k :: Nat) -instance MaximumLikelihood Standard Bernoulli where - mle _ bls = fromList Bernoulli [mean $ toDouble <$> bls] - where toDouble True = 1 - toDouble False = 0 +-- Poisson Distribution -- -instance Legendre Natural Bernoulli where - potential p = log $ 1 + exp (coordinate 0 p) - potentialDifferentials p = fromList (Tangent p) [logistic $ coordinate 0 p] +-- | Returns a sample from a Poisson distribution with the given rate. +samplePoisson :: Double -> Random Int +samplePoisson lmda = Random R.uniform >>= renew 0 + where l = exp (-lmda) + renew k p + | p <= l = return k + | otherwise = do + u <- Random R.uniform + renew (k+1) (p*u) -instance Legendre Mixture Bernoulli where - potential p = - let eta = coordinate 0 p - in logit eta * eta - log (1 / (1 - eta)) - potentialDifferentials p = fromList (Tangent p) [logit $ coordinate 0 p] +-- | The 'Manifold' of 'Poisson' distributions. The 'Source' coordinate is the +-- rate of the Poisson distribution. +data Poisson -instance ExponentialFamily Bernoulli where - baseMeasure _ _ = 1 - sufficientStatistic Bernoulli True = fromList Bernoulli [1] - sufficientStatistic Bernoulli False = fromList Bernoulli [0] +-- von Mises -- -instance Riemannian Natural Bernoulli where - metric p = - let tht = coordinate 0 p - stht = logistic tht - in fromList (Tensor (Tangent p) (Tangent p)) [stht * (1-stht)] +-- | The 'Manifold' of 'VonMises' distributions. The 'Source' coordinates are +-- the mean and concentration. +data VonMises -instance Transition Standard Mixture Bernoulli where - transition = breakChart -instance Transition Mixture Standard Bernoulli where - transition = breakChart +--- Internal --- -instance Transition Standard Natural Bernoulli where - transition = potentialMapping . chart Mixture . transition -instance Transition Natural Standard Bernoulli where - transition = transition . potentialMapping +binomialLogBaseMeasure0 :: (KnownNat n) => Proxy n -> Proxy (Binomial n) -> Int -> Double +binomialLogBaseMeasure0 prxyn _ = logChoose (natValInt prxyn) -instance Generative Natural Bernoulli where - generate = standardGenerate +--- Instances --- --- Binomial Distribution -- -newtype Binomial = Binomial { binomialTrials :: Int } deriving (Eq, Read, Show) +-- Bernoulli Distribution -- -instance Manifold Binomial where - dimension _ = 1 +instance Manifold Bernoulli where + type Dimension Bernoulli = 1 -instance Statistical Binomial where - type SampleSpace Binomial = [Int] - sampleSpace (Binomial n) = [0..n] +instance Statistical Bernoulli where + type (SamplePoint Bernoulli) = Bool -instance Generative Standard Binomial where - generate p = do - let n = binomialTrials $ manifold p - bls <- replicateM n . bernoulli . head $ listCoordinates p - return $ sum [ if bl then 1 else 0 | bl <- bls ] +instance Discrete Bernoulli where + type Cardinality Bernoulli = 2 + sampleSpace _ = [True,False] -instance AbsolutelyContinuous Standard Binomial where - density p k = - let n = binomialTrials $ manifold p - [c] = listCoordinates p - in choose n k * c^k * (1 - c)^(n-k) +instance ExponentialFamily Bernoulli where + logBaseMeasure _ _ = 0 + sufficientStatistic True = Point $ S.singleton 1 + sufficientStatistic False = Point $ S.singleton 0 -instance Legendre Natural Binomial where - potential p = - let n = fromIntegral . binomialTrials $ manifold p - tht = coordinate 0 p - in n * log (1 + exp tht) - potentialDifferentials p = - let n = fromIntegral . binomialTrials $ manifold p - in fromList (Tangent p) [n * logistic (coordinate 0 p)] +type instance PotentialCoordinates Bernoulli = Natural +instance Legendre Bernoulli where + potential p = log $ 1 + exp (S.head $ coordinates p) -instance Legendre Mixture Binomial where - potential p = - let n = fromIntegral . binomialTrials $ manifold p - eta = coordinate 0 p - in eta * log (eta / (n - eta)) - n * log (n / (n - eta)) - potentialDifferentials p = - let n = fromIntegral . binomialTrials $ manifold p - eta = coordinate 0 p - in fromList (Tangent p) [log $ eta / (n - eta) ] +--instance {-# OVERLAPS #-} KnownNat k => Legendre (Replicated k Bernoulli) where +-- potential p = S.sum . S.map (log . (1 +) . exp) $ coordinates p -instance ExponentialFamily Binomial where - baseMeasure (Binomial n) = choose n - sufficientStatistic s k = fromList s [fromIntegral k] +instance Transition Natural Mean Bernoulli where + transition = Point . S.map logistic . coordinates -instance Transition Standard Natural Binomial where - transition = potentialMapping . chart Mixture . transition +instance DuallyFlat Bernoulli where + dualPotential p = + let eta = S.head $ coordinates p + in logit eta * eta - log (1 / (1 - eta)) -instance Transition Natural Standard Binomial where - transition = chart Standard . transition . potentialMapping +instance Transition Mean Natural Bernoulli where + transition = Point . S.map logit . coordinates -instance Transition Standard Mixture Binomial where - transition p = breakChart $ alterCoordinates (* (fromIntegral . binomialTrials $ manifold p)) p +instance Riemannian Natural Bernoulli where + metric p = + let stht = logistic . S.head $ coordinates p + in Point . S.singleton $ stht * (1-stht) + flat p p' = + let stht = logistic . S.head $ coordinates p + in breakPoint $ (stht * (1-stht)) .> p' -instance Transition Mixture Standard Binomial where - transition p = breakChart $ alterCoordinates (/ (fromIntegral . binomialTrials $ manifold p)) p +instance {-# OVERLAPS #-} KnownNat k => Riemannian Natural (Replicated k Bernoulli) where + metric = error "Do not call metric on a replicated manifold" + flat p p' = + let sthts = S.map ((\stht -> stht * (1-stht)) . logistic) $ coordinates p + dp = S.zipWith (*) sthts $ coordinates p' + in Point dp --- Categorical Distribution -- +instance {-# OVERLAPS #-} KnownNat k => Riemannian Mean (Replicated k Bernoulli) where + metric = error "Do not call metric on a replicated manifold" + sharp p dp = + let sthts' = S.map (\stht -> stht * (1-stht)) $ coordinates p + p' = S.zipWith (*) sthts' $ coordinates dp + in Point p' -newtype Categorical s = Categorical s deriving (Show,Eq,Read) --- | A 'Categorical' distribution where the probability of the last category is --- given by the normalization constraint. +instance Transition Source Mean Bernoulli where + transition = breakPoint -generateCategorical :: [k] -> Coordinates -> RandST s k --- | Takes a weighted list of elements representing a probability mass function, and --- returns a sample from the Categorical distribution. -generateCategorical ks0 cs0 = do - c0 <- uniform - return $ findProbability ks0 cs0 c0 - where findProbability ks cs c - | C.null cs = head ks - | c < C.head cs = head ks - | otherwise = findProbability (tail ks) (C.tail cs) (c - C.head cs) +instance Transition Mean Source Bernoulli where + transition = breakPoint -instance Discrete s => Manifold (Categorical s) where - dimension (Categorical s) = length (elements s) - 1 +instance Transition Source Natural Bernoulli where + transition = transition . toMean -instance Discrete s => Statistical (Categorical s) where - type SampleSpace (Categorical s) = s - sampleSpace (Categorical ks) = ks +instance Transition Natural Source Bernoulli where + transition = transition . toMean -instance Discrete s => Generative Standard (Categorical s) where - generate p = generateCategorical (samples $ manifold p) (coordinates p) +instance (Transition c Source Bernoulli) => Generative c Bernoulli where + samplePoint p = Random (R.bernoulli . S.head . coordinates $ toSource p) -instance Discrete s => AbsolutelyContinuous Standard (Categorical s) where - density p k - | idx == dimension (manifold p) = 1 - C.sum cs - | otherwise = cs C.! idx - where cs = coordinates p - idx = fromMaybe (error "attempted to calculate density of non-categorical element") - $ elemIndex k (samples $ manifold p) +instance Transition Mean c Bernoulli => MaximumLikelihood c Bernoulli where + mle = transition . averageSufficientStatistic -instance Discrete s => MaximumLikelihood Standard (Categorical s) where - mle m ks0' = fromIntegral (length ks0') /> fromList m (builder $ samples m) - where builder ks - | null $ tail ks = [] - | otherwise = - let k = head ks - kn = length $ filter (== k) ks0' - in fromIntegral kn : builder (tail ks) +instance LogLikelihood Natural Bernoulli Bool where + logLikelihood = exponentialFamilyLogLikelihood + logLikelihoodDifferential = exponentialFamilyLogLikelihoodDifferential -instance Discrete s => Legendre Natural (Categorical s) where - potential p = log $ 1 + C.sum (exp $ coordinates p) - potentialDifferentials p = - let exps = exp $ coordinates p - nrm = 1 + C.sum exps - in nrm /> fromCoordinates (Tangent p) exps +instance AbsolutelyContinuous Source Bernoulli where + densities sb bs = + let p = S.head $ coordinates sb + in [ if b then p else 1 - p | b <- bs ] -instance Discrete s => Legendre Mixture (Categorical s) where - potential p = - let cs = coordinates p - scs = 1 - C.sum cs - in C.sum (C.zipWith (*) cs $ log cs) + scs * log scs - potentialDifferentials p = - let ps = coordinates p - nrm = 1 - C.sum ps - in fromCoordinates (Tangent p) (log $ C.map (/nrm) ps) +instance AbsolutelyContinuous Mean Bernoulli where + densities = densities . toSource -instance Discrete s => ExponentialFamily (Categorical s) where - baseMeasure _ _ = 1 - sufficientStatistic m k = fromCoordinates m $ C.generate (dimension m) (\j -> if i == j then 1 else 0) - where ks = samples m - i = fromMaybe (error "Categorical distribution given uncategorized element") $ elemIndex k ks +instance AbsolutelyContinuous Natural Bernoulli where + logDensities = exponentialFamilyLogDensities -instance Discrete s => Transition Standard Mixture (Categorical s) where - transition = breakChart +-- Binomial Distribution -- -instance Discrete s => Transition Mixture Standard (Categorical s) where - transition = breakChart +instance KnownNat n => Manifold (Binomial n) where + type Dimension (Binomial n) = 1 -instance Discrete s => Transition Standard Natural (Categorical s) where - transition = potentialMapping . chart Mixture . transition +instance KnownNat n => Statistical (Binomial n) where + type SamplePoint (Binomial n) = Int -instance Discrete s => Transition Natural Standard (Categorical s) where - transition = transition . potentialMapping +instance KnownNat n => Discrete (Binomial n) where + type Cardinality (Binomial n) = n + 1 + sampleSpace prx = [0..binomialSampleSpace prx] --- Curved Categorical Distribution -- +instance KnownNat n => ExponentialFamily (Binomial n) where + logBaseMeasure = binomialLogBaseMeasure0 Proxy + sufficientStatistic = Point . S.singleton . fromIntegral -newtype CurvedCategorical s = CurvedCategorical s deriving (Show,Eq,Read) +type instance PotentialCoordinates (Binomial n) = Natural -instance Discrete s => Manifold (CurvedCategorical s) where - dimension = length . samples +instance KnownNat n => Legendre (Binomial n) where + potential p = + let n = fromIntegral $ binomialTrials p + tht = S.head $ coordinates p + in n * log (1 + exp tht) -instance Discrete s => Statistical (CurvedCategorical s) where - type SampleSpace (CurvedCategorical s) = s - sampleSpace (CurvedCategorical s) = s +instance KnownNat n => Transition Natural Mean (Binomial n) where + transition p = + let n = fromIntegral $ binomialTrials p + in Point . S.singleton $ n * logistic (S.head $ coordinates p) -instance Discrete s => Generative Standard (CurvedCategorical s) where - generate p = generateCategorical (samples $ manifold p) (coordinates p) +instance KnownNat n => DuallyFlat (Binomial n) where + dualPotential p = + let n = fromIntegral $ binomialTrials p + eta = S.head $ coordinates p + in eta * log (eta / (n - eta)) - n * log (n / (n - eta)) -instance Discrete s => AbsolutelyContinuous Standard (CurvedCategorical s) where - density p k = cs C.! idx - where ks = samples $ manifold p - cs = coordinates p - idx = fromMaybe (error "attempted to calculate density of non-categorical element") - $ elemIndex k ks +instance KnownNat n => Transition Mean Natural (Binomial n) where + transition p = + let n = fromIntegral $ binomialTrials p + eta = S.head $ coordinates p + in Point . S.singleton . log $ eta / (n - eta) --- Poisson Distribution -- +instance KnownNat n => Transition Source Natural (Binomial n) where + transition = transition . toMean -generatePoisson :: Double -> RandST s Int --- | Returns a sample from a Poisson distribution with the given rate. -generatePoisson rt = - uniform >>= renew 0 - where l = exp (-rt) - renew k p - | p <= l = return k - | otherwise = do - u <- uniform - renew (k+1) (p*u) +instance KnownNat n => Transition Natural Source (Binomial n) where + transition = transition . toMean -data Poisson = Poisson deriving (Eq, Read, Show) +instance KnownNat n => Transition Source Mean (Binomial n) where + transition p = + let n = fromIntegral $ binomialTrials p + in breakPoint $ n .> p -instance Manifold Poisson where - dimension _ = 1 +instance KnownNat n => Transition Mean Source (Binomial n) where + transition p = + let n = fromIntegral $ binomialTrials p + in breakPoint $ n /> p -instance Statistical Poisson where - type SampleSpace Poisson = NaturalNumbers - sampleSpace _ = NaturalNumbers +instance (KnownNat n, Transition c Source (Binomial n)) => Generative c (Binomial n) where + samplePoint p0 = do + let p = toSource p0 + n = binomialTrials p + rb = Random (R.bernoulli . S.head $ coordinates p) + bls <- replicateM n rb + return $ sum [ if bl then 1 else 0 | bl <- bls ] -instance Generative Standard Poisson where - generate d = generatePoisson . C.head $ coordinates d +instance KnownNat n => AbsolutelyContinuous Source (Binomial n) where + densities p ks = + let n = binomialTrials p + c = S.head $ coordinates p + in [ choose n k * c^k * (1 - c)^(n-k) | k <- ks ] -instance AbsolutelyContinuous Standard Poisson where - density d k = - let ps = coordinates d - lmda = C.head ps - in lmda^k / factorial k * exp (-lmda) +instance KnownNat n => AbsolutelyContinuous Mean (Binomial n) where + densities = densities . toSource -instance MaximumLikelihood Standard Poisson where - mle _ xs = fromList Poisson . (:[]) . mean $ fromIntegral <$> xs +instance KnownNat n => AbsolutelyContinuous Natural (Binomial n) where + logDensities = exponentialFamilyLogDensities -instance ExponentialFamily Poisson where - sufficientStatistic Poisson = fromCoordinates Poisson . C.singleton . fromIntegral - baseMeasure _ k = recip $ factorial k +instance (KnownNat n, Transition Mean c (Binomial n)) => MaximumLikelihood c (Binomial n) where + mle = transition . averageSufficientStatistic -instance Legendre Natural Poisson where - potential p = exp $ coordinate 0 p - potentialDifferentials p = fromCoordinates (Tangent p) . exp $ coordinates p +instance KnownNat n => LogLikelihood Natural (Binomial n) Int where + logLikelihood = exponentialFamilyLogLikelihood + logLikelihoodDifferential = exponentialFamilyLogLikelihoodDifferential -instance Legendre Mixture Poisson where - potential p = - let eta = coordinate 0 p - in eta * log eta - eta - potentialDifferentials p = fromCoordinates (Tangent p) . log $ coordinates p -instance Riemannian Natural Poisson where - metric p = - let tht = coordinate 0 p - in fromList (Tensor (Tangent p) (Tangent p)) [exp tht] +-- Categorical Distribution -- -instance Transition Standard Natural Poisson where - transition = transition . chart Mixture . transition +instance KnownNat n => Manifold (Categorical n) where + type Dimension (Categorical n) = n -instance Transition Natural Standard Poisson where - transition = transition . potentialMapping +instance KnownNat n => Statistical (Categorical n) where + type SamplePoint (Categorical n) = Int -instance Transition Standard Mixture Poisson where - transition = breakChart +instance KnownNat n => Discrete (Categorical n) where + type Cardinality (Categorical n) = n + sampleSpace prx = [0..dimension prx] -instance Transition Mixture Standard Poisson where - transition = breakChart +instance KnownNat n => ExponentialFamily (Categorical n) where + logBaseMeasure _ _ = 0 + sufficientStatistic e = Point $ S.generate (\i -> if finiteInt i == (fromEnum e-1) then 1 else 0) -instance Generative Natural Poisson where - generate = standardGenerate +type instance (PotentialCoordinates (Categorical n)) = Natural --- Normal Distribution -- +instance KnownNat n => Legendre (Categorical n) where + --potential (Point cs) = log $ 1 + S.sum (S.map exp cs) + potential = logSumExp . B.cons 0 . boxCoordinates -data Normal = Normal deriving (Show,Eq,Read) +instance KnownNat n => Transition Natural Mean (Categorical n) where + transition p = + let exps = S.map exp $ coordinates p + nrm = 1 + S.sum exps + in nrm /> Point exps -instance Manifold Normal where - dimension _ = 2 +instance KnownNat n => DuallyFlat (Categorical n) where + dualPotential (Point cs) = + let sc = 1 - S.sum cs + in S.sum (S.map entropyFun cs) + entropyFun sc + where entropyFun 0 = 0 + entropyFun x = x * log x -instance Statistical Normal where - type SampleSpace Normal = Continuum - sampleSpace _ = Continuum +instance KnownNat n => Transition Mean Natural (Categorical n) where + transition (Point xs) = + let nrm = 1 - S.sum xs + in Point . log $ S.map (/nrm) xs -instance Generative Standard Normal where - generate p = - let [mu,vr] = listCoordinates p - in normal mu $ sqrt vr +instance Transition Source Mean (Categorical n) where + transition = breakPoint -instance AbsolutelyContinuous Standard Normal where - density p x = - let [mu,vr] = listCoordinates p - in recip (sqrt $ vr*2*pi) * exp (negate $ (x - mu) ** 2 / (2*vr)) +instance Transition Mean Source (Categorical n) where + transition = breakPoint -instance MaximumLikelihood Standard Normal where - mle _ xs = - let (mu,vr) = meanVariance $ C.fromList xs - in fromList Normal [mu,vr] +instance KnownNat n => Transition Source Natural (Categorical n) where + transition = transition . toMean -instance ExponentialFamily Normal where - sufficientStatistic Normal x = fromList Normal [x,x**2] - baseMeasure _ _ = recip . sqrt $ 2 * pi +instance KnownNat n => Transition Natural Source (Categorical n) where + transition = transition . toMean -instance Legendre Natural Normal where - potential p = - let [tht0,tht1] = listCoordinates p - in -(tht0^2 / (4*tht1)) - 0.5 * log(-2*tht1) - potentialDifferentials p = - let [tht0,tht1] = listCoordinates p - dv = tht0/tht1 - in fromList (Tangent p) [-0.5*dv, 0.25 * dv^2 - 0.5/tht1] +instance (KnownNat n, Transition c Source (Categorical n)) + => Generative c (Categorical n) where + samplePoint p0 = + let p = toSource p0 + in sampleCategorical $ coordinates p -instance Legendre Mixture Normal where - potential p = - let [eta0,eta1] = listCoordinates p - in -0.5 * log(eta1 - eta0^2) - 1/2 - potentialDifferentials p = - let [eta0,eta1] = listCoordinates p - dff = eta0^2 - eta1 - in fromList (Tangent p) [-eta0 / dff, 0.5 / dff] +instance (KnownNat n, Transition Mean c (Categorical n)) + => MaximumLikelihood c (Categorical n) where + mle = transition . averageSufficientStatistic -instance Riemannian Natural Normal where - metric p = - let [tht1,tht2] = listCoordinates p - in fromList (Tensor (Tangent p) (Tangent p)) - [-1/(2*tht2),tht1/(2*tht2^2),tht1/(2*tht2^2),(-tht1^2 + tht2)/(2*tht2^3) ] +instance KnownNat n => LogLikelihood Natural (Categorical n) Int where + logLikelihood = exponentialFamilyLogLikelihood + logLikelihoodDifferential = exponentialFamilyLogLikelihoodDifferential -instance Riemannian Standard Normal where - metric p = - let [_,vr] = listCoordinates p - in fromList (Tensor (Tangent p) (Tangent p)) [recip vr,0,0,recip $ 2*vr^2] -instance Transition Standard Mixture Normal where - transition p = - let [mu,vr] = listCoordinates p - in fromList Normal [mu, vr + mu^2] +instance KnownNat n => AbsolutelyContinuous Source (Categorical n) where + densities (Point ps) es = do + e <- es + let ek = fromEnum e + p0 = 1 - S.sum ps + return $ if ek == 0 + then p0 + else S.unsafeIndex ps $ ek - 1 -instance Transition Mixture Standard Normal where - transition p = - let [eta0,eta1] = listCoordinates p - in fromList Normal [eta0, eta1 - eta0^2] +instance KnownNat n => AbsolutelyContinuous Mean (Categorical n) where + densities = densities . toSource -instance Transition Standard Natural Normal where - transition p = - let [mu,vr] = listCoordinates p - in fromList Normal [mu / vr, negate . recip $ 2 * vr] +instance KnownNat n => AbsolutelyContinuous Natural (Categorical n) where + logDensities = exponentialFamilyLogDensities -instance Transition Natural Standard Normal where - transition p = - let [tht0,tht1] = listCoordinates p - in fromList Normal [-0.5 * tht0 / tht1, negate . recip $ 2 * tht1] +-- Dirichlet Distribution -- -instance Generative Natural Normal where - generate = standardGenerate +instance KnownNat k => Manifold (Dirichlet k) where + type Dimension (Dirichlet k) = k --- MeanNormal Distribution -- +instance KnownNat k => Statistical (Dirichlet k) where + type SamplePoint (Dirichlet k) = S.Vector k Double -data MeanNormal = MeanNormal Double deriving (Show,Eq,Read) +instance (KnownNat k, Transition c Source (Dirichlet k)) + => Generative c (Dirichlet k) where + samplePoint p0 = do + let alphs = boxCoordinates $ toSource p0 + G.convert <$> Random (R.dirichlet alphs) -instance Manifold MeanNormal where - dimension _ = 1 +instance KnownNat k => ExponentialFamily (Dirichlet k) where + logBaseMeasure _ = negate . S.sum + sufficientStatistic xs = Point $ S.map log xs +logMultiBeta :: KnownNat k => S.Vector k Double -> Double +logMultiBeta alphs = + S.sum (S.map GSL.lngamma alphs) - GSL.lngamma (S.sum alphs) -instance Statistical MeanNormal where - type SampleSpace MeanNormal = Continuum - sampleSpace _ = Continuum +logMultiBetaDifferential :: KnownNat k => S.Vector k Double -> S.Vector k Double +logMultiBetaDifferential alphs = + S.map (subtract (GSL.psi $ S.sum alphs) . GSL.psi) alphs -instance Generative Standard MeanNormal where - generate p = do - let (MeanNormal vr) = manifold p - normal (coordinate 0 p) $ sqrt vr +type instance PotentialCoordinates (Dirichlet k) = Natural -instance AbsolutelyContinuous Standard MeanNormal where - density p = - let (MeanNormal vr) = manifold p - mu = coordinate 0 p - in density . chart Standard $ fromList Normal [mu,vr] +instance KnownNat k => Legendre (Dirichlet k) where + potential = logMultiBeta . coordinates -instance MaximumLikelihood Standard MeanNormal where - mle mnrm xs = fromList mnrm [mean xs] +instance KnownNat k => Transition Natural Mean (Dirichlet k) where + transition = Point . logMultiBetaDifferential . coordinates -instance Legendre Natural MeanNormal where - potential p = - let (MeanNormal vr) = manifold p - in 0.5 * vr * coordinate 0 p^2 - potentialDifferentials p = - let (MeanNormal vr) = manifold p - in fromList (Tangent p) [vr * coordinate 0 p] +instance KnownNat k => AbsolutelyContinuous Source (Dirichlet k) where + densities p xss = do + xs <- xss + let alphs = coordinates p + prds = S.product $ S.zipWith (**) xs $ S.map (subtract 1) alphs + return $ prds / exp (logMultiBeta alphs) -instance Legendre Mixture MeanNormal where - potential p = - let (MeanNormal vr) = manifold p - in 0.5 / vr * coordinate 0 p^2 - potentialDifferentials p = - let (MeanNormal vr) = manifold p - in fromList (Tangent p) [coordinate 0 p / vr] +instance KnownNat k => AbsolutelyContinuous Natural (Dirichlet k) where + logDensities = exponentialFamilyLogDensities -instance ExponentialFamily MeanNormal where - sufficientStatistic mnrm x = fromList mnrm [x] - baseMeasure (MeanNormal vr) x = (exp . negate $ 0.5 * x^2 / vr) / sqrt (2*pi*vr) +instance KnownNat k => LogLikelihood Natural (Dirichlet k) (S.Vector k Double) where + logLikelihood = exponentialFamilyLogLikelihood + logLikelihoodDifferential = exponentialFamilyLogLikelihoodDifferential -instance Riemannian Natural MeanNormal where - metric p = - let (MeanNormal vr) = manifold p - in fromList (Tensor (Tangent p) (Tangent p)) [vr] +instance KnownNat k => Transition Source Natural (Dirichlet k) where + transition = breakPoint -instance Transition Standard Natural MeanNormal where - transition = potentialMapping . chart Mixture . breakChart +instance KnownNat k => Transition Natural Source (Dirichlet k) where + transition = breakPoint -instance Transition Natural Standard MeanNormal where - transition = breakChart . potentialMapping +-- Poisson Distribution -- -instance Transition Standard Mixture MeanNormal where - transition = breakChart +instance Manifold Poisson where + type Dimension Poisson = 1 -instance Transition Mixture Standard MeanNormal where - transition = breakChart +instance Statistical Poisson where + type SamplePoint Poisson = Int --- Multivariate Normal -- +instance ExponentialFamily Poisson where + sufficientStatistic = Point . S.singleton . fromIntegral + logBaseMeasure _ k = negate $ logFactorial k -data MultivariateNormal = MultivariateNormal { sampleSpaceDimension :: Int } deriving (Eq, Read, Show) +type instance PotentialCoordinates Poisson = Natural -generateMultivariateNormal :: C.Vector Double -> M.Matrix Double -> RandST s (C.Vector Double) --- | Samples from a multivariate Normal. -generateMultivariateNormal mus rtsgma = do - nrms <- C.replicateM n $ normal 0 1 - return $ mus + (M.#>) rtsgma nrms - where n = C.length mus +instance Legendre Poisson where + potential = exp . S.head . coordinates -muSigmaToMultivariateNormal :: C.Vector Double -> M.Matrix Double -> Standard :#: MultivariateNormal --- | Generates a multivariateNormal by way of a covariance matrix i.e. by taking --- the square root. -muSigmaToMultivariateNormal mus sgma = - fromCoordinates (MultivariateNormal $ C.length mus) $ mus C.++ M.flatten sgma +instance Transition Natural Mean Poisson where + transition = Point . exp . coordinates -splitCoordinates :: c :#: MultivariateNormal -> (Coordinates, M.Matrix Double) -splitCoordinates p = - let (MultivariateNormal n) = manifold p - (mus,sgms) = C.splitAt n $ coordinates p - in (mus,M.reshape n sgms) +instance DuallyFlat Poisson where + dualPotential (Point xs) = + let eta = S.head xs + in eta * log eta - eta -instance Manifold MultivariateNormal where - dimension (MultivariateNormal n) = n + n^2 +instance Transition Mean Natural Poisson where + transition = Point . log . coordinates -instance Statistical MultivariateNormal where - type SampleSpace MultivariateNormal = Euclidean - sampleSpace (MultivariateNormal n) = Euclidean n +instance Transition Source Natural Poisson where + transition = transition . toMean -instance Generative Standard MultivariateNormal where - generate p = - let n = sampleSpaceDimension $ manifold p - (mus,sds) = C.splitAt n $ coordinates p - in generateMultivariateNormal mus $ M.reshape n sds +instance Transition Natural Source Poisson where + transition = transition . toMean -instance AbsolutelyContinuous Standard MultivariateNormal where - density p xs = - let n = sampleSpaceDimension $ manifold p - (mus,sgma) = splitCoordinates p - flx = M.sqrtm sgma - in recip ((2*pi)**(fromIntegral n / 2) * M.det flx) - * exp (-0.5 * ((M.tr (M.inv sgma) M.#> C.zipWith (-) xs mus) `M.dot` C.zipWith (-) xs mus)) +instance Transition Source Mean Poisson where + transition = breakPoint -instance MaximumLikelihood Standard MultivariateNormal where - mle _ xss = - let n = fromIntegral $ length xss - mus = recip (fromIntegral n) * sum xss - sgma = recip (fromIntegral $ n - 1) - * sum (map (\xs -> let xs' = xs - mus in M.outer xs' xs') xss) - in muSigmaToMultivariateNormal mus sgma +instance Transition Mean Source Poisson where + transition = breakPoint -instance ExponentialFamily MultivariateNormal where - sufficientStatistic m x = fromCoordinates m $ x C.++ M.flatten (M.outer x x) - baseMeasure (MultivariateNormal n) _ = (2*pi)**(-fromIntegral n/2) +instance (Transition c Source Poisson) => Generative c Poisson where + samplePoint = samplePoisson . S.head . coordinates . toSource -instance Legendre Natural MultivariateNormal where - potential p = - let (tmu,tsgma) = splitCoordinates p - invtsgma = M.inv tsgma - in -0.25 * M.dot tmu (invtsgma M.#> tmu) - 0.5 * log(M.det $ M.scale (-2) tsgma) - potentialDifferentials p = - let (tmu,tsgma) = splitCoordinates p - invtsgma = M.inv tsgma - invapp = M.app invtsgma tmu - in fromCoordinates (Tangent p) $ (-0.5 * invapp) - C.++ M.flatten (M.scale (-0.5) invtsgma + M.scale 0.25 (M.outer invapp invapp)) +instance AbsolutelyContinuous Source Poisson where + densities (Point xs) ks = do + k <- ks + let lmda = S.head xs + return $ lmda^k / factorial k * exp (-lmda) -instance Legendre Mixture MultivariateNormal where - potential p = - let (mmu,msgma) = splitCoordinates p - --n = fromIntegral . sampleSpaceDimension $ manifold p - in -0.5 * (1 + M.dot mmu (M.inv msgma M.#> mmu)) - 0.5 * log (M.det msgma) - potentialDifferentials p = - let (mmu,msgma) = splitCoordinates p - invmsgma' = M.inv $ M.outer mmu mmu - msgma - in fromCoordinates (Tangent p) $ (negate invmsgma' M.#> mmu) C.++ M.flatten (M.scale 0.5 invmsgma') +instance AbsolutelyContinuous Mean Poisson where + densities = densities . toSource -instance Transition Standard Natural MultivariateNormal where - transition p = - let (mu,sgma) = splitCoordinates p - invsgma = M.inv sgma - in fromCoordinates (manifold p) $ (invsgma M.#> mu) C.++ M.flatten (M.scale (-0.5) invsgma) +instance AbsolutelyContinuous Natural Poisson where + logDensities = exponentialFamilyLogDensities -instance Transition Natural Standard MultivariateNormal where - transition p = - let (emu,esgma) = splitCoordinates p - invesgma = M.inv esgma - in fromCoordinates (manifold p) $ M.scale 0.5 (invesgma M.#> emu) C.++ M.flatten (M.scale 0.5 invesgma) +instance Transition Mean c Poisson => MaximumLikelihood c Poisson where + mle = transition . averageSufficientStatistic -instance Transition Standard Mixture MultivariateNormal where - transition p = - let (mu,sgma) = splitCoordinates p - in fromCoordinates (manifold p) $ mu C.++ M.flatten (sgma + M.outer mu mu) +instance LogLikelihood Natural Poisson Int where + logLikelihood = exponentialFamilyLogLikelihood + logLikelihoodDifferential = exponentialFamilyLogLikelihoodDifferential -instance Transition Mixture Standard MultivariateNormal where - transition p = - let (mmu,msgma) = splitCoordinates p - in fromCoordinates (manifold p) $ mmu C.++ M.flatten (msgma -M.outer mmu mmu) +-- VonMises -- +instance Manifold VonMises where + type Dimension VonMises = 2 -{- ---- Graveyard --- +instance Statistical VonMises where + type SamplePoint VonMises = Double +instance Generative Source VonMises where + samplePoint p@(Point cs) = do + let (mu,kap0) = S.toPair cs + kap = max kap0 1e-5 + tau = 1 + sqrt (1 + 4 * square kap) + rho = (tau - sqrt (2*tau))/(2*kap) + r = (1 + square rho) / (2 * rho) + u1 <- Random R.uniform + u2 <- Random R.uniform + u3 <- Random R.uniform + let z = cos (pi * u1) + f = (1 + r * z)/(r + z) + c = kap * (r - f) + if log (c / u2) + 1 - c < 0 + then samplePoint p + else return . toPi $ signum (u3 - 0.5) * acos f + mu -functionToCategorical :: Double -> Double -> Int -> (Double -> Double) -> Standard :#: Categorical Double --- | Takes range information in the form of a minimum, maximum, and sample count, --- and a function which represents an unnomralized pdf, and returns a normalized list of --- pairs (x,f(x)) over the specified range such that the sum of the f(x)s is 1. --- --- In principle, f should be strictly positive, but this is not checked. -functionToCategorical mn mx n f = - let (ks,fks) = unzip $ discretizeFunction mn mx n f - in recip (sum fks) .> fromList (Categorical ks) fks +instance AbsolutelyContinuous Source VonMises where + densities p xs = do + let (mu,kp) = S.toPair $ coordinates p + x <- xs + return $ exp (kp * cos (x - mu)) / (2*pi * GSL.bessel_I0 kp) --- Exponential Distribution -- +instance LogLikelihood Natural VonMises Double where + logLikelihood = exponentialFamilyLogLikelihood + logLikelihoodDifferential = exponentialFamilyLogLikelihoodDifferential -data Exponential = Exponential deriving (Eq,Read,Show) +type instance PotentialCoordinates VonMises = Natural -instance Manifold Exponential where - dimension _ = 1 +instance Legendre VonMises where + potential p = + let kp = snd . S.toPair . coordinates $ toSource p + in log $ GSL.bessel_I0 kp -type instance SampleSpace Exponential = Continuum +instance Transition Natural Mean VonMises where + transition p = + let kp = snd . S.toPair . coordinates $ toSource p + in breakPoint $ (GSL.bessel_I1 kp / (GSL.bessel_I0 kp * kp)) .> p -instance Statistical Exponential where - sampleSpace _ = Continuum +instance AbsolutelyContinuous Natural VonMises where + logDensities = exponentialFamilyLogDensities -instance Generative Standard Exponential where - generate = exponential . C.head . coordinates +instance Generative Natural VonMises where + samplePoint = samplePoint . toSource -instance AbsolutelyContinuous Standard Exponential where - density p x = - let lmda = C.head $ coordinates p - in lmda * exp (negate $ lmda * x) +instance ExponentialFamily VonMises where + sufficientStatistic tht = Point $ S.doubleton (cos tht) (sin tht) + logBaseMeasure _ _ = -log(2 * pi) -instance MaximumLikelihood Standard Exponential where - mle _ xs = chart Standard . fromList Exponential . (:[]) . recip . mean $ xs +instance Transition Source Natural VonMises where + transition (Point cs) = + let (mu,kap) = S.toPair cs + in Point $ S.doubleton (kap * cos mu) (kap * sin mu) -instance Legendre Natural Exponential where - potential p = negate . log . negate $ coordinate 0 p - potentialDifferentials p = fromCoordinates (Tangent p) . negate $ coordinates p +instance Transition Natural Source VonMises where + transition (Point cs) = + let (tht0,tht1) = S.toPair cs + in Point $ S.doubleton (toPi $ atan2 tht1 tht0) (sqrt $ square tht0 + square tht1) -instance Legendre Mixture Exponential where - potential p = 1 - log eta - potentialDifferentials p = +instance Transition Source Mean VonMises where + transition = toMean . toNatural -instance ExponentialFamily Exponential where - sufficientStatistic Exponential = fromCoordinates Exponential . C.singleton - baseMeasure _ _ = 1 -instance Transition Standard Natural Exponential where - transition = breakChart . alterCoordinates negate +--- Location Shape --- -instance Transition Natural Standard Exponential where - transition = breakChart . alterCoordinates negate +instance (Statistical l, Manifold s) => Statistical (LocationShape l s) where + type SamplePoint (LocationShape l s) = SamplePoint l --} +instance (Manifold l, Manifold s) => Translation (LocationShape l s) l where + (>+>) yz y' = + let (y,z) = split yz + in join (y + y') z + anchor = fst . split + +type instance PotentialCoordinates (LocationShape l s) = Natural + +instance ( Statistical l, Statistical s , Product (LocationShape l s) + , Storable (SamplePoint s), SamplePoint l ~ SamplePoint s + , AbsolutelyContinuous c (LocationShape l s), KnownNat n) + => AbsolutelyContinuous c (LocationShape (Replicated n l) (Replicated n s)) where + logDensities lss xs = + let (l,s) = split lss + ls = splitReplicated l + ss = splitReplicated s + lss' :: c # Replicated n (LocationShape l s) + lss' = joinReplicated $ S.zipWith join ls ss + in logDensities lss' xs + + +instance (KnownNat n, Manifold l, Manifold s) + => Translation (Replicated n (LocationShape l s)) (Replicated n l) where + {-# INLINE (>+>) #-} + (>+>) w z = + let ws = splitReplicated w + zs = splitReplicated z + in joinReplicated $ S.zipWith (>+>) ws zs + {-# INLINE anchor #-} + anchor = mapReplicatedPoint anchor diff --git a/Goal/Probability/Distributions/CoMPoisson.hs b/Goal/Probability/Distributions/CoMPoisson.hs new file mode 100644 --- /dev/null +++ b/Goal/Probability/Distributions/CoMPoisson.hs @@ -0,0 +1,195 @@ +{-# OPTIONS_GHC -fplugin=GHC.TypeLits.KnownNat.Solver -fplugin=GHC.TypeLits.Normalise -fconstraint-solver-iterations=10 #-} +{-# LANGUAGE UndecidableInstances #-} + +-- | Implementation of Conway-Maxwell Poisson distributions (CoMPoisson). +-- () CoMPoisson distributions generalize Poisson distributions with +-- a shape parameter that can concentrate or disperse the underlying Poisson +-- distribution. +module Goal.Probability.Distributions.CoMPoisson + ( + -- * CoMPoisson + CoMPoisson + , CoMShape + -- ** Numerics + , comPoissonLogPartitionSum + , comPoissonExpectations + ) where + +-- Package -- + +import Goal.Core +import Goal.Geometry + +import Goal.Probability.Statistical +import Goal.Probability.ExponentialFamily +import Goal.Probability.Distributions + +import qualified Goal.Core.Vector.Storable as S +import qualified System.Random.MWC as R + + +--- Analysis --- + +--- CoMPoisson Distribution --- + +-- | A type for storing the shape of a 'CoMPoisson' distribution. +data CoMShape + +-- | The 'Manifold' of 'CoMPoisson' distributions. The 'Source' coordinates of the +-- 'CoMPoisson' are the mode $\mu$ and the "pseudo-precision" parameter $\nu$, such that $\mu / \nu$ is approximately the variance of the distribution. +type CoMPoisson = LocationShape Poisson CoMShape + +-- | Approximates the log-partition function of the given CoMPoisson +-- distribution up to the specified precision. +comPoissonLogPartitionSum :: Double -> Natural # CoMPoisson -> Double +{-# INLINE comPoissonLogPartitionSum #-} +comPoissonLogPartitionSum eps np = + let (tht1,tht2) = S.toPair $ coordinates np + in fst $ comPoissonLogPartitionSum0 eps tht1 tht2 + +-- | Approximates the expectations of functions given the natural parameters of +-- a CoM-Poisson distribution. +comPoissonExpectations + :: KnownNat n + => Double + -> (Int -> S.Vector n Double) + -> Natural # CoMPoisson + -> S.Vector n Double +{-# INLINE comPoissonExpectations #-} +comPoissonExpectations eps f np = + let (tht1,tht2) = S.toPair $ coordinates np + (lgprt,ln) = comPoissonLogPartitionSum0 eps tht1 tht2 + js = [0..ln] + dns = exp . subtract lgprt <$> unnormalizedLogDensities np js + in sum $ zipWith S.scale dns (f <$> js) + +-- | Approximates the mean mparameters of a CoM-Poisson distribution. +comPoissonMeans :: Double -> Natural # CoMPoisson -> Mean # CoMPoisson +{-# INLINE comPoissonMeans #-} +comPoissonMeans eps cp = + let ss :: Int -> Mean # CoMPoisson + ss = sufficientStatistic + in Point $ comPoissonExpectations eps (coordinates . ss) cp + + +--- Internal --- + + +comPoissonSequence :: Double -> Double -> [Double] +comPoissonSequence tht1 tht2 = + [ tht1 * fromIntegral j + logFactorial j *tht2 | (j :: Int) <- [0..] ] + +comPoissonLogPartitionSum0 :: Double -> Double -> Double -> (Double, Int) +{-# INLINE comPoissonLogPartitionSum0 #-} +comPoissonLogPartitionSum0 eps tht1 tht2 = + let md = floor $ comPoissonSmoothMode tht1 tht2 + (hdsqs,tlsqs) = splitAt md $ comPoissonSequence tht1 tht2 + mx = tht1 * fromIntegral md + logFactorial md *tht2 + ehdsqs = exp . subtract mx <$> hdsqs + etlsqs = exp . subtract mx <$> tlsqs + sqs' = ehdsqs ++ takeWhile (> eps) etlsqs + in ((+ mx) . log1p . subtract 1 $ sum sqs' , length sqs') + +comPoissonSmoothMode :: Double -> Double -> Double +comPoissonSmoothMode tht1 tht2 = exp (tht1/negate tht2) + +--comPoissonApproximateMean :: Double -> Double -> Double +--comPoissonApproximateMean mu nu = +-- mu + 1/(2*nu) - 0.5 +-- +--comPoissonApproximateVariance :: Double -> Double -> Double +--comPoissonApproximateVariance mu nu = mu / nu + +overDispersedEnvelope :: Double -> Double -> Double -> Double +overDispersedEnvelope p mu nu = + let mnm1 = 1 - p + flrd = max 0 . floor $ mu / (mnm1**recip nu) + nmr = mu**(nu * fromIntegral flrd) + dmr = (mnm1^flrd) * (factorial flrd ** nu) + in recip p * nmr / dmr + +underDispersedEnvelope :: Double -> Double -> Double +underDispersedEnvelope mu nu = + let fmu = floor mu + in (mu ^ fmu / factorial fmu)** (nu - 1) + +sampleOverDispersed :: Double -> Double -> Double -> Double -> Random Int +sampleOverDispersed p bnd0 mu nu = do + u0 <- Random R.uniform + let y' = max 0 . floor $ logBase (1 - p) u0 + nmr = (mu^y' / factorial y')**nu + dmr = bnd0 * (1-p)^y' * p + alph = nmr/dmr + u <- Random R.uniform + if isNaN alph + then error "NaN in sampling CoMPoisson: Parameters out of bounds" + else if u <= alph + then return y' + else sampleOverDispersed p bnd0 mu nu + +sampleUnderDispersed :: Double -> Double -> Double -> Random Int +sampleUnderDispersed bnd0 mu nu = do + let psn :: Source # Poisson + psn = Point $ S.singleton mu + y' <- samplePoint psn + let alph0 = mu^y' / factorial y' + alph = alph0**nu / (bnd0*alph0) + u <- Random R.uniform + if u <= alph + then return y' + else sampleUnderDispersed bnd0 mu nu + +sampleCoMPoisson :: Int -> Double -> Double -> Random [Int] +sampleCoMPoisson n mu nu + | nu >= 1 = + let bnd0 = underDispersedEnvelope mu nu + in replicateM n $ sampleUnderDispersed bnd0 mu nu + | otherwise = + let p = 2*nu / (2*mu*nu + 1 + nu) + bnd0 = overDispersedEnvelope p mu nu + in replicateM n $ sampleOverDispersed p bnd0 mu nu + + +-- Instances -- + + +instance ExponentialFamily CoMPoisson where + sufficientStatistic k = fromTuple (fromIntegral k, logFactorial k) + logBaseMeasure _ _ = 0 + +type instance PotentialCoordinates CoMPoisson = Natural + +instance Legendre CoMPoisson where + potential = + comPoissonLogPartitionSum 1e-16 + +instance AbsolutelyContinuous Natural CoMPoisson where + logDensities = exponentialFamilyLogDensities + +instance Transition Source Natural CoMPoisson where + transition p = + let (mu,nu) = S.toPair $ coordinates p + in fromTuple (nu * log mu, -nu) + +instance Transition Natural Source CoMPoisson where + transition p = + let (tht1,tht2) = S.toPair $ coordinates p + in fromTuple (exp (-tht1/tht2), -tht2) + +instance (Transition c Source CoMPoisson) => Generative c CoMPoisson where + sample n p = do + let (mu,nu) = S.toPair . coordinates $ toSource p + in sampleCoMPoisson n mu nu + +instance Transition Natural Mean CoMPoisson where + transition = comPoissonMeans 1e-16 + +instance Transition Source Mean CoMPoisson where + transition = toMean . toNatural + +instance LogLikelihood Natural CoMPoisson Int where + logLikelihood = exponentialFamilyLogLikelihood + logLikelihoodDifferential = exponentialFamilyLogLikelihoodDifferential + +instance Manifold CoMShape where + type Dimension CoMShape = 1 diff --git a/Goal/Probability/Distributions/Gaussian.hs b/Goal/Probability/Distributions/Gaussian.hs new file mode 100644 --- /dev/null +++ b/Goal/Probability/Distributions/Gaussian.hs @@ -0,0 +1,560 @@ +{-# OPTIONS_GHC -fplugin=GHC.TypeLits.KnownNat.Solver -fplugin=GHC.TypeLits.Normalise -fconstraint-solver-iterations=10 #-} +{-# LANGUAGE UndecidableInstances,TypeApplications #-} + +-- | Various instances of statistical manifolds, with a focus on exponential +-- families. In the documentation we use \(X\) to indicate a random variable +-- with the distribution being documented. +module Goal.Probability.Distributions.Gaussian + ( -- * Univariate + Normal + , NormalMean + , NormalVariance + -- * Multivariate + , MVNMean + , MVNCovariance + , MultivariateNormal + , multivariateNormalCorrelations + , bivariateNormalConfidenceEllipse + , splitMultivariateNormal + , splitMeanMultivariateNormal + , splitNaturalMultivariateNormal + , joinMultivariateNormal + , joinMeanMultivariateNormal + , joinNaturalMultivariateNormal + -- * Linear Models + , SimpleLinearModel + , LinearModel + ) where + +-- Package -- + +import Goal.Core +import Goal.Probability.Statistical +import Goal.Probability.ExponentialFamily +import Goal.Probability.Distributions + +import Goal.Geometry + +import qualified Goal.Core.Vector.Storable as S +import qualified Goal.Core.Vector.Generic as G + +import qualified System.Random.MWC.Distributions as R + +-- Normal Distribution -- + +-- | The Mean of a normal distribution. When used as a distribution itself, it +-- is a Normal distribution with unit variance. +data NormalMean + +-- | The variance of a normal distribution. +data NormalVariance + +-- | The 'Manifold' of 'Normal' distributions. The 'Source' coordinates are the +-- mean and the variance. +type Normal = LocationShape NormalMean NormalVariance + +-- | The Mean of a normal distribution. When used as a distribution itself, it +-- is a Normal distribution with unit variance. +data MVNMean (n :: Nat) + +-- | The variance of a normal distribution. +data MVNCovariance (n :: Nat) + +-- | Linear models are linear functions with additive Guassian noise. +type LinearModel n k = Affine Tensor (MVNMean n) (MultivariateNormal n) (MVNMean k) + +-- | Linear models are linear functions with additive Guassian noise. +type SimpleLinearModel = Affine Tensor NormalMean Normal NormalMean + +-- Multivariate Normal -- + +-- | The 'Manifold' of 'MultivariateNormal' distributions. The 'Source' +-- coordinates are the (vector) mean and the covariance matrix. For the +-- coordinates of a multivariate normal distribution, the elements of the mean +-- come first, and then the elements of the covariance matrix in row major +-- order. +-- +-- Note that we only store the lower triangular elements of the covariance +-- matrix, to better reflect the true dimension of a MultivariateNormal +-- Manifold. In short, be careful when using 'join' and 'split' to access the +-- values of the Covariance matrix, and consider using the specific instances +-- for MVNs. +type MultivariateNormal (n :: Nat) = LocationShape (MVNMean n) (MVNCovariance n) + +-- | Split a MultivariateNormal into its Means and Covariance matrix. +splitMultivariateNormal + :: KnownNat n + => Source # MultivariateNormal n + -> (S.Vector n Double, S.Matrix n n Double) +splitMultivariateNormal mvn = + let (mu,cvr) = split mvn + in (coordinates mu, S.fromLowerTriangular $ coordinates cvr) + +-- | Join a covariance matrix into a MultivariateNormal. +joinMultivariateNormal + :: KnownNat n + => S.Vector n Double + -> S.Matrix n n Double + -> Source # MultivariateNormal n +joinMultivariateNormal mus sgma = + join (Point mus) (Point $ S.lowerTriangular sgma) + +-- | Split a MultivariateNormal into its Means and Covariance matrix. +splitMeanMultivariateNormal + :: KnownNat n + => Mean # MultivariateNormal n + -> (S.Vector n Double, S.Matrix n n Double) +splitMeanMultivariateNormal mvn = + let (mu,cvr) = split mvn + in (coordinates mu, S.fromLowerTriangular $ coordinates cvr) + +-- | Join a covariance matrix into a MultivariateNormal. +joinMeanMultivariateNormal + :: KnownNat n + => S.Vector n Double + -> S.Matrix n n Double + -> Mean # MultivariateNormal n +joinMeanMultivariateNormal mus sgma = + join (Point mus) (Point $ S.lowerTriangular sgma) + +-- | Split a MultivariateNormal into the precision weighted means and (-0.5*) +-- Precision matrix. Note that this performs an easy to miss computation for +-- converting the natural parameters in our reduced representation of MVNs into +-- the full precision matrix. +splitNaturalMultivariateNormal + :: KnownNat n + => Natural # MultivariateNormal n + -> (S.Vector n Double, S.Matrix n n Double) +splitNaturalMultivariateNormal np = + let (nmu,cvrs) = split np + nmu0 = coordinates nmu + nsgma0' = (/2) . S.fromLowerTriangular $ coordinates cvrs + nsgma0 = nsgma0' + S.diagonalMatrix (S.takeDiagonal nsgma0') + in (nmu0, nsgma0) + +-- | Joins a MultivariateNormal out of the precision weighted means and (-0.5) +-- Precision matrix. Note that this performs an easy to miss computation for +-- converting the full precision Matrix into the reduced, EF representation we use here. +joinNaturalMultivariateNormal + :: KnownNat n + => S.Vector n Double + -> S.Matrix n n Double + -> Natural # MultivariateNormal n +joinNaturalMultivariateNormal nmu0 nsgma0 = + let nmu = Point nmu0 + diag = S.diagonalMatrix $ S.takeDiagonal nsgma0 + in join nmu . Point . S.lowerTriangular $ 2*nsgma0 - diag + +-- | Confidence elipses for bivariate normal distributions. +bivariateNormalConfidenceEllipse + :: Int + -> Double + -> Source # MultivariateNormal 2 + -> [(Double,Double)] +bivariateNormalConfidenceEllipse nstps prcnt nrm = + let (mu,cvr) = splitMultivariateNormal nrm + chl = S.withMatrix (S.scale prcnt) $ S.unsafeCholesky cvr + xs = range 0 (2*pi) nstps + sxs = [ S.fromTuple (cos x, sin x) | x <- xs ] + in S.toPair . (mu +) <$> S.matrixMap chl sxs + +-- | Computes the correlation matrix of a 'MultivariateNormal' distribution. +multivariateNormalCorrelations + :: KnownNat k + => Source # MultivariateNormal k + -> S.Matrix k k Double +multivariateNormalCorrelations mnrm = + let cvrs = snd $ splitMultivariateNormal mnrm + sds = S.map sqrt $ S.takeDiagonal cvrs + sdmtx = S.outerProduct sds sds + in G.Matrix $ S.zipWith (/) (G.toVector cvrs) (G.toVector sdmtx) + +multivariateNormalLogBaseMeasure + :: forall n . (KnownNat n) + => Proxy (MultivariateNormal n) + -> S.Vector n Double + -> Double +multivariateNormalLogBaseMeasure _ _ = + let n = natValInt (Proxy :: Proxy n) + in -fromIntegral n/2 * log (2*pi) + +mvnMeanLogBaseMeasure + :: forall n . (KnownNat n) + => Proxy (MVNMean n) + -> S.Vector n Double + -> Double +mvnMeanLogBaseMeasure _ x = + let n = natValInt (Proxy :: Proxy n) + in -fromIntegral n/2 * log pi - S.dotProduct x x / 2 + +-- | samples a multivariateNormal by way of a covariance matrix i.e. by taking +-- the square root. +sampleMultivariateNormal + :: KnownNat n + => Source # MultivariateNormal n + -> Random (S.Vector n Double) +sampleMultivariateNormal p = do + let (mus,sgma) = splitMultivariateNormal p + nrms <- S.replicateM $ Random (R.normal 0 1) + let rtsgma = S.matrixRoot sgma + return $ mus + S.matrixVectorMultiply rtsgma nrms + + +--- Internal --- + + +--- Instances --- + + +-- NormalMean Distribution -- + +instance Manifold NormalMean where + type Dimension NormalMean = 1 + +instance Statistical NormalMean where + type SamplePoint NormalMean = Double + +instance ExponentialFamily NormalMean where + sufficientStatistic x = singleton x + logBaseMeasure _ x = -square x/2 - sqrt (2*pi) + +type instance PotentialCoordinates NormalMean = Natural + +instance Transition Mean Natural NormalMean where + transition = breakPoint + +instance Transition Mean Source NormalMean where + transition = breakPoint + +instance Transition Source Natural NormalMean where + transition = breakPoint + +instance Transition Source Mean NormalMean where + transition = breakPoint + +instance Transition Natural Mean NormalMean where + transition = breakPoint + +instance Transition Natural Source NormalMean where + transition = breakPoint + +instance Legendre NormalMean where + potential (Point cs) = + let tht = S.head cs + in square tht / 2 + +instance LogLikelihood Natural NormalMean Double where + logLikelihood = exponentialFamilyLogLikelihood + logLikelihoodDifferential = exponentialFamilyLogLikelihoodDifferential + + +-- Normal Shape -- + + +instance Manifold NormalVariance where + type Dimension NormalVariance = 1 + + +-- Normal Distribution -- + +instance ExponentialFamily Normal where + sufficientStatistic x = + Point . S.doubleton x $ x**2 + logBaseMeasure _ _ = -1/2 * log (2 * pi) + +type instance PotentialCoordinates Normal = Natural + +instance Legendre Normal where + potential (Point cs) = + let (tht0,tht1) = S.toPair cs + in -(square tht0 / (4*tht1)) - 0.5 * log(-2*tht1) + +instance Transition Natural Mean Normal where + transition p = + let (tht0,tht1) = S.toPair $ coordinates p + dv = tht0/tht1 + in Point $ S.doubleton (-0.5*dv) (0.25 * square dv - 0.5/tht1) + +instance DuallyFlat Normal where + dualPotential (Point cs) = + let (eta0,eta1) = S.toPair cs + in -0.5 * log(eta1 - square eta0) - 1/2 + +instance Transition Mean Natural Normal where + transition p = + let (eta0,eta1) = S.toPair $ coordinates p + dff = eta1 - square eta0 + in Point $ S.doubleton (eta0 / dff) (-0.5 / dff) + +instance Riemannian Natural Normal where + metric p = + let (tht0,tht1) = S.toPair $ coordinates p + d00 = -1/(2*tht1) + d01 = tht0/(2*square tht1) + d11 = 0.5*(1/square tht1 - square tht0 / (tht1^(3 :: Int))) + in Point $ S.doubleton d00 d01 S.++ S.doubleton d01 d11 + +instance Riemannian Mean Normal where + metric p = + let (eta0,eta1) = S.toPair $ coordinates p + eta02 = square eta0 + dff2 = square $ eta1 - eta02 + d00 = (dff2 + 2 * eta02) / dff2 + d01 = -eta0 / dff2 + d11 = 0.5 / dff2 + in Point $ S.doubleton d00 d01 S.++ S.doubleton d01 d11 + +-- instance Riemannian Source Normal where +-- metric p = +-- let (_,vr) = S.toPair $ coordinates p +-- in Point $ S.doubleton (recip vr) 0 S.++ S.doubleton 0 (recip $ 2*square vr) + +instance Transition Source Mean Normal where + transition (Point cs) = + let (mu,vr) = S.toPair cs + in Point . S.doubleton mu $ vr + square mu + +instance Transition Mean Source Normal where + transition (Point cs) = + let (eta0,eta1) = S.toPair cs + in Point . S.doubleton eta0 $ eta1 - square eta0 + +instance Transition Source Natural Normal where + transition (Point cs) = + let (mu,vr) = S.toPair cs + in Point $ S.doubleton (mu / vr) (negate . recip $ 2 * vr) + +instance Transition Natural Source Normal where + transition (Point cs) = + let (tht0,tht1) = S.toPair cs + in Point $ S.doubleton (-0.5 * tht0 / tht1) (negate . recip $ 2 * tht1) + +instance (Transition c Source Normal) => Generative c Normal where + samplePoint p = + let (Point cs) = toSource p + (mu,vr) = S.toPair cs + in Random $ R.normal mu (sqrt vr) + +instance AbsolutelyContinuous Source Normal where + densities (Point cs) xs = do + let (mu,vr) = S.toPair cs + x <- xs + return $ recip (sqrt $ vr*2*pi) * exp (negate $ (x - mu) ** 2 / (2*vr)) + +instance AbsolutelyContinuous Mean Normal where + densities = densities . toSource + +instance AbsolutelyContinuous Natural Normal where + logDensities = exponentialFamilyLogDensities + +instance Transition Mean c Normal => MaximumLikelihood c Normal where + mle = transition . averageSufficientStatistic + +instance LogLikelihood Natural Normal Double where + logLikelihood = exponentialFamilyLogLikelihood + logLikelihoodDifferential = exponentialFamilyLogLikelihoodDifferential + + +-- MVNMean -- + +instance KnownNat n => Manifold (MVNMean n) where + type Dimension (MVNMean n) = n + +instance (KnownNat n) => Statistical (MVNMean n) where + type SamplePoint (MVNMean n) = S.Vector n Double + +instance KnownNat n => ExponentialFamily (MVNMean n) where + sufficientStatistic x = Point x + logBaseMeasure = mvnMeanLogBaseMeasure + +type instance PotentialCoordinates (MVNMean n) = Natural + +-- MVNCovariance -- + +instance (KnownNat n, KnownNat (Triangular n)) => Manifold (MVNCovariance n) where + type Dimension (MVNCovariance n) = Triangular n + +-- Multivariate Normal -- + +instance (KnownNat n, KnownNat (Triangular n)) + => AbsolutelyContinuous Source (MultivariateNormal n) where + densities mvn xs = do + let (mu,sgma) = splitMultivariateNormal mvn + n = fromIntegral $ natValInt (Proxy @ n) + scl = (2*pi)**(-n/2) * S.determinant sgma**(-1/2) + isgma = S.pseudoInverse sgma + x <- xs + let dff = x - mu + expval = S.dotProduct dff $ S.matrixVectorMultiply isgma dff + return $ scl * exp (-expval / 2) + +instance (KnownNat n, KnownNat (Triangular n), Transition c Source (MultivariateNormal n)) + => Generative c (MultivariateNormal n) where + samplePoint = sampleMultivariateNormal . toSource + +instance KnownNat n => Transition Source Natural (MultivariateNormal n) where + transition p = + let (mu,sgma) = splitMultivariateNormal p + invsgma = S.pseudoInverse sgma + in joinNaturalMultivariateNormal (S.matrixVectorMultiply invsgma mu) $ (-0.5) * invsgma + +instance KnownNat n => Transition Natural Source (MultivariateNormal n) where + transition p = + let (nmu,nsgma) = splitNaturalMultivariateNormal p + insgma = (-0.5) * S.pseudoInverse nsgma + in joinMultivariateNormal (S.matrixVectorMultiply insgma nmu) insgma + +instance KnownNat n => LogLikelihood Natural (MultivariateNormal n) (S.Vector n Double) where + logLikelihood = exponentialFamilyLogLikelihood + logLikelihoodDifferential = exponentialFamilyLogLikelihoodDifferential + + +instance (KnownNat n, KnownNat (Triangular n)) => ExponentialFamily (MultivariateNormal n) where + sufficientStatistic xs = Point $ xs S.++ S.lowerTriangular (S.outerProduct xs xs) + averageSufficientStatistic xs = Point $ average xs S.++ S.lowerTriangular ( S.averageOuterProduct $ zip xs xs ) + logBaseMeasure = multivariateNormalLogBaseMeasure + +type instance PotentialCoordinates (MultivariateNormal n) = Natural + +instance (KnownNat n, KnownNat (Triangular n)) => Legendre (MultivariateNormal n) where + potential p = + let (nmu,nsgma) = splitNaturalMultivariateNormal p + insgma = S.pseudoInverse nsgma + in -0.25 * S.dotProduct nmu (S.matrixVectorMultiply insgma nmu) + -0.5 * (log . S.determinant . negate $ 2 * nsgma) + +instance (KnownNat n, KnownNat (Triangular n)) => Transition Natural Mean (MultivariateNormal n) where + transition = toMean . toSource + +instance (KnownNat n, KnownNat (Triangular n)) => DuallyFlat (MultivariateNormal n) where + dualPotential p = + let sgma = snd . splitMultivariateNormal $ toSource p + n = natValInt (Proxy @ n) + lndet = fromIntegral n*log (2*pi*exp 1) + log (S.determinant sgma) + in -0.5 * lndet + +instance (KnownNat n, KnownNat (Triangular n)) => Transition Mean Natural (MultivariateNormal n) where + transition = toNatural . toSource + +instance KnownNat n => Transition Source Mean (MultivariateNormal n) where + transition p = + let (mu,sgma) = splitMultivariateNormal p + in joinMeanMultivariateNormal mu $ sgma + S.outerProduct mu mu + +instance KnownNat n => Transition Mean Source (MultivariateNormal n) where + transition p = + let (mu,scnds) = splitMeanMultivariateNormal p + in joinMultivariateNormal mu $ scnds - S.outerProduct mu mu + +instance (KnownNat n, KnownNat (Triangular n)) => AbsolutelyContinuous Natural (MultivariateNormal n) where + logDensities = exponentialFamilyLogDensities + +instance (KnownNat n, Transition Mean c (MultivariateNormal n)) + => MaximumLikelihood c (MultivariateNormal n) where + mle = transition . averageSufficientStatistic + +--instance KnownNat n => MaximumLikelihood Source (MultivariateNormal n) where +-- mle _ xss = +-- let n = fromIntegral $ length xss +-- mus = recip (fromIntegral n) * sum xss +-- sgma = recip (fromIntegral $ n - 1) +-- * sum (map (\xs -> let xs' = xs - mus in M.outer xs' xs') xss) +-- in joinMultivariateNormal mus sgma + +-- Linear Models + +instance ( KnownNat n, KnownNat k) + => Transition Natural Source (Affine Tensor (MVNMean n) (MultivariateNormal n) (MVNMean k)) where + transition nfa = + let (mvn,nmtx) = split nfa + (nmu,nsg) = splitNaturalMultivariateNormal mvn + invsg = -2 * nsg + ssg = S.inverse invsg + smu = S.matrixVectorMultiply ssg nmu + smvn = joinMultivariateNormal smu ssg + smtx = S.matrixMatrixMultiply ssg $ toMatrix nmtx + in join smvn $ fromMatrix smtx + +instance ( KnownNat n, KnownNat k) + => Transition Source Natural (Affine Tensor (MVNMean n) (MultivariateNormal n) (MVNMean k)) where + transition lmdl = + let (smvn,smtx) = split lmdl + (smu,ssg) = splitMultivariateNormal smvn + invsg = S.inverse ssg + nmu = S.matrixVectorMultiply invsg smu + nsg = -0.5 * invsg + nmtx = S.matrixMatrixMultiply invsg $ toMatrix smtx + nmvn = joinNaturalMultivariateNormal nmu nsg + in join nmvn $ fromMatrix nmtx + +instance ( KnownNat n, KnownNat k) + => Transition Natural Source (Affine Tensor (MVNMean n) (Replicated n Normal) (MVNMean k)) where + transition nfa = + let (nnrms,nmtx) = split nfa + (nmu,nsg) = splitReplicatedProduct nnrms + nmvn = joinNaturalMultivariateNormal (coordinates nmu) $ S.diagonalMatrix (coordinates nsg) + nlm :: Natural # LinearModel n k + nlm = join nmvn nmtx + (smvn,smtx) = split $ transition nlm + (smu,ssg) = splitMultivariateNormal smvn + snrms = joinReplicatedProduct (Point smu) (Point $ S.takeDiagonal ssg) + in join snrms smtx + +instance ( KnownNat n, KnownNat k) + => Transition Source Natural (Affine Tensor (MVNMean n) (Replicated n Normal) (MVNMean k)) where + transition sfa = + let (snrms,smtx) = split sfa + (smu,ssg) = S.toPair . S.toColumns . S.fromRows . S.map coordinates $ splitReplicated snrms + smvn = joinMultivariateNormal smu $ S.diagonalMatrix ssg + slm :: Source # LinearModel n k + slm = join smvn smtx + (nmvn,nmtx) = split $ transition slm + (nmu,nsg) = splitNaturalMultivariateNormal nmvn + nnrms = joinReplicated $ S.zipWith (curry fromTuple) nmu $ S.takeDiagonal nsg + in join nnrms nmtx + +instance Transition Natural Source (Affine Tensor NormalMean Normal NormalMean) where + transition nfa = + let nfa' :: Natural # LinearModel 1 1 + nfa' = breakPoint nfa + sfa' :: Source # LinearModel 1 1 + sfa' = transition nfa' + in breakPoint sfa' + +instance Transition Source Natural (Affine Tensor NormalMean Normal NormalMean) where + transition sfa = + let sfa' :: Source # LinearModel 1 1 + sfa' = breakPoint sfa + nfa' :: Natural # LinearModel 1 1 + nfa' = transition sfa' + in breakPoint nfa' + + + +--instance ( KnownNat n, KnownNat k) +-- => Transition Natural Source (Affine Tensor (MVNMean n) (Replicated n Normal) (MVNMean k)) where +-- transition nfa = +-- let (nnrms,nmtx) = split nfa +-- (nmu,nsg) = S.toPair . S.toColumns . S.fromRows . S.map coordinates +-- $ splitReplicated nnrms +-- invsg = -2 * nsg +-- ssg = recip invsg +-- smu = nmu / invsg +-- snrms = joinReplicated $ S.zipWith (curry fromTuple) smu ssg +-- smtx = S.matrixMatrixMultiply (S.diagonalMatrix ssg) $ toMatrix nmtx +-- in join snrms $ fromMatrix smtx + +--instance ( KnownNat n, KnownNat k) +-- => Transition Source Natural (Affine Tensor (MVNMean n) (Replicated n Normal) (MVNMean k)) where +-- transition sfa = +-- let (snrms,smtx) = split sfa +-- (smu,ssg) = S.toPair . S.toColumns . S.fromRows . S.map coordinates +-- $ splitReplicated snrms +-- invsg = recip ssg +-- nmu = invsg * smu +-- nsg = -0.5 * invsg +-- nmtx = S.matrixMatrixMultiply (S.diagonalMatrix invsg) $ toMatrix smtx +-- nnrms = joinReplicated $ S.zipWith (curry fromTuple) nmu nsg +-- in join nnrms $ fromMatrix nmtx + + diff --git a/Goal/Probability/ExponentialFamily.hs b/Goal/Probability/ExponentialFamily.hs --- a/Goal/Probability/ExponentialFamily.hs +++ b/Goal/Probability/ExponentialFamily.hs @@ -1,13 +1,30 @@ -module Goal.Probability.ExponentialFamily ( - -- * Exponential Families - ExponentialFamily (sufficientStatistic, baseMeasure) - , sufficientStatisticN - -- ** Dual Parameters - , Natural (Natural) - , Mixture (Mixture) - -- ** Divergence - , klDivergence +{-# LANGUAGE UndecidableInstances,TypeApplications #-} +-- | Definitions for working with exponential families. +module Goal.Probability.ExponentialFamily + ( -- * Exponential Families + ExponentialFamily (sufficientStatistic, averageSufficientStatistic, logBaseMeasure) + , LegendreExponentialFamily + , DuallyFlatExponentialFamily + , exponentialFamilyLogDensities + , unnormalizedLogDensities + -- ** Coordinate Systems + , Natural + , Mean + , Source + -- ** Coordinate Transforms + , toNatural + , toMean + , toSource + -- ** Entropies , relativeEntropy + , crossEntropy + -- ** Differentials + , relativeEntropyDifferential + , stochasticRelativeEntropyDifferential + , stochasticInformationProjectionDifferential + -- *** Maximum Likelihood Instances + , exponentialFamilyLogLikelihood + , exponentialFamilyLogLikelihoodDifferential ) where --- Imports --- @@ -17,82 +34,198 @@ import Goal.Probability.Statistical +import Goal.Core import Goal.Geometry +import qualified Goal.Core.Vector.Storable as S +import Foreign.Storable + --- Exponential Families --- --- | A 'Statistical' 'Manifold' is a member of the 'ExponentialFamily' if we can --- specify a 'sufficientStatistic' of fixed length. Defining the 'baseMeasure' --- is also necessary in order to render unique the 'Natural' and 'Mixture' --- parameterizations. +-- | A parameterization which represents the standard or typical parameterization of +-- the given manifold, e.g. the Poisson rate or Normal mean and standard deviation. +data Source + +-- | A parameterization in terms of the natural parameters of an exponential family. +data Natural + +-- | A parameterization in terms of the mean 'sufficientStatistic' of an exponential family. +data Mean + +instance Primal Natural where + type Dual Natural = Mean + +instance Primal Mean where + type Dual Mean = Natural + +-- | Expresses an exponential family distribution in 'Natural' coordinates. +toNatural :: (Transition c Natural x) => c # x -> Natural # x +toNatural = transition + +-- | Expresses an exponential family distribution in 'Mean' coordinates. +toMean :: (Transition c Mean x) => c # x -> Mean # x +toMean = transition + +-- | Expresses an exponential family distribution in 'Source' coordinates. +toSource :: (Transition c Source x) => c # x -> Source # x +toSource = transition + +-- | An 'ExponentialFamily' is a 'Statistical' 'Manifold' \( \mathcal M \) +-- determined by a fixed-length 'sufficientStatistic' \(s_i\) and a +-- 'logBaseMeasure' \(\mu\). Each distribution \(P \in \mathcal M\) may then be +-- identified with 'Natural' parameters \(\theta_i\) such that +-- \(p(x) \propto e^{\sum_{i=1}^n \theta_i s_i(x)}\mu(x)\). 'ExponentialFamily' +-- distributions theoretically have a 'Riemannian' geometry, with 'metric' +-- 'Tensor' given by the Fisher information metric. However, not all +-- distributions (e.g. the von Mises distribution) afford closed-form +-- expressions for all the relevant structures. +class Statistical x => ExponentialFamily x where + sufficientStatistic :: SamplePoint x -> Mean # x + averageSufficientStatistic :: Sample x -> Mean # x + averageSufficientStatistic = average . map sufficientStatistic + logBaseMeasure :: Proxy x -> SamplePoint x -> Double + +-- | When the log-partition function and its derivative of the given +-- 'ExponentialFamily' may be computed in closed-form, then we refer to it as a +-- 'LegendreExponentialFamily'. -- --- 'ExponentialFamily' distributions theoretically have a 'Riemannian' geometry --- given by the Fisher information metric, given rise to the 'DualChart' system --- of 'Natural' and 'Mixture'. A 'Point' on the 'ExponentialFamily' 'Manifold' in --- one of these dual coordinates is assumed to be equipped the corresponding --- dual connection. Under this assumption, we take the 'Manifold' itself to be --- self-dual to simplify types. -class (Statistical m, Legendre Natural m, Legendre Mixture m) => ExponentialFamily m where - sufficientStatistic :: m -> Sample m -> Mixture :#: m - baseMeasure :: m -> Sample m -> Double +-- Note that the log-partition function is the 'potential' of the 'Legendre' +-- class, and its derivative maps 'Natural' coordinates to 'Mean' coordinates. +type LegendreExponentialFamily x = + ( PotentialCoordinates x ~ Natural, Legendre x, ExponentialFamily x + , Transition (PotentialCoordinates x) (Dual (PotentialCoordinates x)) x ) -sufficientStatisticN :: ExponentialFamily m => m -> [Sample m] -> Mixture :#: m --- | The sufficient statistic of N iid random variables. -sufficientStatisticN m xs = - fromIntegral (length xs) /> foldr1 (<+>) (sufficientStatistic m <$> xs) +-- | When additionally, the (negative) entropy and its derivative of the given +-- 'ExponentialFamily' may be computed in closed-form, then we refer to it as a +-- 'DuallyFlatExponentialFamily'. +-- +-- Note that the negative entropy is the 'dualPotential' of the 'DuallyFlat' class, +-- and its derivative maps 'Mean' coordinates to 'Natural' coordinates. +type DuallyFlatExponentialFamily x = + ( LegendreExponentialFamily x, DuallyFlat x + , Transition (Dual (PotentialCoordinates x)) (PotentialCoordinates x) x ) -klDivergence - :: (ExponentialFamily m, Transition c Natural m, Transition d Mixture m) - => c :#: m -> d :#: m -> Double -klDivergence q p = divergence (chart Natural $ transition q) (chart Mixture $ transition p) +-- | The relative entropy \(D(P \parallel Q)\), also known as the KL-divergence. +-- This is simply the 'canonicalDivergence' with its arguments flipped. +relativeEntropy :: DuallyFlatExponentialFamily x => Mean # x -> Natural # x -> Double +relativeEntropy = flip canonicalDivergence -relativeEntropy - :: (ExponentialFamily m, Transition c Mixture m, Transition d Natural m) - => c :#: m -> d :#: m -> Double -relativeEntropy p q = klDivergence q p +-- | A function for computing the cross-entropy, which is the relative entropy +-- plus the entropy of the first distribution. +crossEntropy :: DuallyFlatExponentialFamily x => Mean # x -> Natural # x -> + Double +crossEntropy mp nq = potential nq - (mp <.> nq) --- | A parameterization in terms of the natural coordinates of an exponential family. -data Natural = Natural +-- | The differential of the relative entropy with respect to the 'Natural' parameters of +-- the second argument. +relativeEntropyDifferential :: LegendreExponentialFamily x => Mean # x -> Natural # x -> Mean # x +relativeEntropyDifferential mp nq = transition nq - mp --- | A representation in terms of the mean sufficient statistics of an exponential family. -data Mixture = Mixture +-- | Monte Carlo estimate of the differential of the relative entropy with +-- respect to the 'Natural' parameters of the second argument, based on samples from +-- the two distributions. +stochasticRelativeEntropyDifferential + :: ExponentialFamily x + => Sample x -- ^ True Samples + -> Sample x -- ^ Model Samples + -> Mean # x -- ^ Differential Estimate +stochasticRelativeEntropyDifferential pxs qxs = + averageSufficientStatistic qxs - averageSufficientStatistic pxs -instance Primal Natural where - type Dual Natural = Mixture +-- | Estimate of the differential of relative entropy with respect to the +-- 'Natural' parameters of the first argument, based a 'Sample' from the first +-- argument and the unnormalized log-density of the second. +stochasticInformationProjectionDifferential + :: ExponentialFamily x + => Natural # x -- ^ Model Distribution + -> Sample x -- ^ Model Samples + -> (SamplePoint x -> Double) -- ^ Unnormalized log-density of target distribution + -> Mean # x -- ^ Differential Estimate +stochasticInformationProjectionDifferential px xs f = + let mxs = sufficientStatistic <$> xs + mys = (\x -> sufficientStatistic x <.> px - f x) <$> xs + ln = fromIntegral $ length xs + mxht = ln /> sum mxs + myht = sum mys / ln + in (ln - 1) /> sum [ (my - myht) .> (mx - mxht) | (mx,my) <- zip mxs mys ] -instance Primal Mixture where - type Dual Mixture = Natural +-- | The density of an exponential family distribution that has an exact +-- expression for the log-partition function. +exponentialFamilyLogDensities + :: (ExponentialFamily x, Legendre x, PotentialCoordinates x ~ Natural) => Natural # x -> Sample x -> [Double] +exponentialFamilyLogDensities p xs = subtract (potential p) <$> unnormalizedLogDensities p xs +-- | The unnormalized log-density of an arbitrary exponential family distribution. +unnormalizedLogDensities :: forall x . ExponentialFamily x => Natural # x -> Sample x -> [Double] +unnormalizedLogDensities p xs = + zipWith (+) (dotMap p $ sufficientStatistic <$> xs) (logBaseMeasure (Proxy @ x) <$> xs) ---- Instances --- +-- | 'logLikelihood' for a 'LegendreExponentialFamily'. +exponentialFamilyLogLikelihood + :: forall x . LegendreExponentialFamily x + => Sample x -> Natural # x -> Double +exponentialFamilyLogLikelihood xs nq = + let mp = averageSufficientStatistic xs + bm = average $ logBaseMeasure (Proxy :: Proxy x) <$> xs + in -potential nq + (mp <.> nq) + bm +-- | 'logLikelihoodDifferential' for a 'LegendreExponentialFamily'. +exponentialFamilyLogLikelihoodDifferential + :: LegendreExponentialFamily x + => Sample x -> Natural # x -> Mean # x +exponentialFamilyLogLikelihoodDifferential xs nq = + let mp = averageSufficientStatistic xs + in mp - transition nq --- Generic -- -instance ExponentialFamily m => MaximumLikelihood Mixture m where - mle = sufficientStatisticN +--- Internal --- -instance ExponentialFamily m => MaximumLikelihood Natural m where - mle m xs = potentialMapping $ sufficientStatisticN m xs +replicatedlogBaseMeasure0 :: (ExponentialFamily x, Storable (SamplePoint x), KnownNat k) + => Proxy x -> Proxy (Replicated k x) -> S.Vector k (SamplePoint x) -> Double +replicatedlogBaseMeasure0 prxym _ xs = S.sum $ S.map (logBaseMeasure prxym) xs + +pairlogBaseMeasure + :: (ExponentialFamily x, ExponentialFamily y) + => Proxy x + -> Proxy y + -> Proxy (x,y) + -> SamplePoint (x,y) + -> Double +pairlogBaseMeasure prxym prxyn _ (xm,xn) = + logBaseMeasure prxym xm + logBaseMeasure prxyn xn + + +--- Instances --- + + -- Replicated -- -instance ExponentialFamily m => ExponentialFamily (Replicated m) where - sufficientStatistic (Replicated m _) xs = - joinReplicated $ sufficientStatistic m <$> xs - baseMeasure (Replicated m _) xs = product $ baseMeasure m <$> xs +instance Transition Natural Natural x where + transition = id --- Fisher Manifolds -- +instance Transition Mean Mean x where + transition = id -instance ExponentialFamily m => AbsolutelyContinuous Natural m where - density p x = - let s = manifold p - in exp ((p <.> sufficientStatistic s x) - potential p) * baseMeasure s x +instance Transition Source Source x where + transition = id -instance ExponentialFamily m => Transition Mixture Natural m where - transition = potentialMapping +instance (ExponentialFamily x, Storable (SamplePoint x), KnownNat k) + => ExponentialFamily (Replicated k x) where + sufficientStatistic xs = joinReplicated $ S.map sufficientStatistic xs + logBaseMeasure = replicatedlogBaseMeasure0 Proxy -instance ExponentialFamily m => Transition Natural Mixture m where - transition = potentialMapping +-- Sum -- + +instance (ExponentialFamily x, ExponentialFamily y) => ExponentialFamily (x,y) where + sufficientStatistic (xm,xn) = + join (sufficientStatistic xm) (sufficientStatistic xn) + logBaseMeasure = pairlogBaseMeasure Proxy Proxy + + +-- Source Coordinates -- + +instance Primal Source where + type Dual Source = Source diff --git a/Goal/Probability/Graphical.hs b/Goal/Probability/Graphical.hs deleted file mode 100644 --- a/Goal/Probability/Graphical.hs +++ /dev/null @@ -1,9 +0,0 @@ -module Goal.Probability.Graphical where - -import Goal.Geometry -import Goal.Probability.ExponentialFamily - --- | A 'Function' from the 'Mixture' 'Coordinates' of one 'ExponentialFamily' to --- another. Fundamental to neural networks of various kinds. -type NaturalFunction = Function Mixture Natural - diff --git a/Goal/Probability/Graphical/Harmonium.hs b/Goal/Probability/Graphical/Harmonium.hs deleted file mode 100644 --- a/Goal/Probability/Graphical/Harmonium.hs +++ /dev/null @@ -1,214 +0,0 @@ --- | Exponential Family 'Harmonium's and gibbs sampling. -module Goal.Probability.Graphical.Harmonium - ( -- * Harmoniums - Harmonium (Harmonium) - -- ** Type Synonyms - , NaturalFunction - -- ** Structural Manipulation - , splitHarmonium - , joinHarmonium - , harmoniumTranspose - -- ** Conditional Distribution Functions - , conditionalLatentDistribution - , conditionalObservableDistribution - , conditionalLatentDistributions - , conditionalObservableDistributions - -- ** Gibbs Sampling - , bulkGibbsSampling - , bulkGibbsSampling0 - -- * Transducers - , buildNormalTransducer - , buildReplicatedNormalTransducer - , modulateTransducerGain - , modulateHarmoniumBelief - ) where - - ---- Imports --- - - --- Goal -- - -import Goal.Geometry - -import Goal.Probability.Statistical -import Goal.Probability.ExponentialFamily -import Goal.Probability.Distributions -import Goal.Probability.Graphical - -import System.Random.MWC.Monad -import qualified Data.Vector.Storable as C - - ---- Types --- - --- | A quadratic function in the product space of two exponential families. -data Harmonium m n = Harmonium m n deriving (Eq, Read, Show) - --- Datatype manipulation -- - -splitHarmonium :: (Manifold m, Manifold n) - => Function c d :#: Harmonium m n -> (d :#: m, Function c d :#: Tensor m n, Dual c :#: n) --- | Splits a 'Harmonium' into its components parts of a 'Tensor' and a pair of biases. -splitHarmonium qdc = - let (Harmonium m n) = manifold qdc - tns = Tensor m n - (mcs,css') = C.splitAt (dimension m) $ coordinates qdc - (mtxcs,ncs) = C.splitAt (dimension tns) css' - in (fromCoordinates m mcs, fromCoordinates tns mtxcs, fromCoordinates n ncs) - -joinHarmonium - :: (Manifold m, Manifold n) => d :#: m -> Function c d :#: Tensor m n -> Dual c :#: n -> Function c d :#: Harmonium m n --- | Assembles a 'Harmonium' out of the components of the quadratic function. -joinHarmonium dm mtx cn = - let (Tensor m n) = manifold mtx - in fromCoordinates (Harmonium m n) $ coordinates dm C.++ coordinates mtx C.++ coordinates cn - -harmoniumTranspose :: (Manifold n, Manifold m, Primal c, Primal d) - => Function c d :#: Harmonium m n -> Function (Dual d) (Dual c) :#: Harmonium n m --- | Transposes the 'Tensor' in the 'Harmonium' and swaps the biases. -harmoniumTranspose qdc = - let (dm,mtx,dn) = splitHarmonium qdc - in joinHarmonium dn (matrixTranspose mtx) dm - - ---- Functions --- - - -conditionalLatentDistributions :: (Manifold m, ExponentialFamily n) - => NaturalFunction :#: Harmonium m n -> [Sample n] -> [Natural :#: m] --- | Calculates the latent distributions given some observations. -conditionalLatentDistributions p os = - let (Harmonium _ n) = manifold p - in p >$> (sufficientStatistic n <$> os) - -conditionalObservableDistributions :: (ExponentialFamily m, Manifold n) - => NaturalFunction :#: Harmonium m n -> [Sample m] -> [Natural :#: n] --- | Calculates the observable distributions given some latent states. -conditionalObservableDistributions p ls = - let (Harmonium m _) = manifold p - in harmoniumTranspose p >$> (sufficientStatistic m <$> ls) - -conditionalLatentDistribution :: (Manifold m, ExponentialFamily n) - => NaturalFunction :#: Harmonium m n -> Sample n -> Natural :#: m --- | Calculates the latent distributions given an observation. -conditionalLatentDistribution p o = - let (Harmonium _ n) = manifold p - in p >.> sufficientStatistic n o - -conditionalObservableDistribution :: (ExponentialFamily m, Manifold n) - => NaturalFunction :#: Harmonium m n -> Sample m -> Natural :#: n --- | Calculates the observable distributions given a latent state. -conditionalObservableDistribution p l = - let (Harmonium m _) = manifold p - in harmoniumTranspose p >.> sufficientStatistic m l - -bulkGibbsSampling - :: (ExponentialFamily m, Generative Natural m, ExponentialFamily n, Generative Natural n) - => Int -> NaturalFunction :#: Harmonium m n -> [Sample n] -> RandST s [[(Sample m, Sample n)]] --- | Returns a Markov chain over the latent and observable states generated by Gibbs sampling. -bulkGibbsSampling k0 p o0s = do - l0s <- mapM generate $ conditionalLatentDistributions p o0s - gbs <- gibbsSampler k0 l0s [] - return $ zip l0s o0s : gbs - where (Harmonium m n) = manifold p - gibbsSampler 0 _ acc = return $ reverse acc - gibbsSampler k ls acc = do - let mls = sufficientStatistic m <$> ls - os' <- mapM generate $ harmoniumTranspose p >$> mls - let mos' = sufficientStatistic n <$> os' - ls' <- mapM generate $ p >$> mos' - gibbsSampler (k-1) ls' (zip ls' os':acc) - -bulkGibbsSampling0 - :: (ExponentialFamily m, Generative Natural m, ExponentialFamily n, Generative Natural n) - => Int -> NaturalFunction :#: Harmonium m n -> [Mixture :#: n] -> RandST s [[(Mixture :#: m, Mixture :#: n)]] --- | Returns a Markov chain over the latent and observable expoential families generated by Gibbs sampling. -bulkGibbsSampling0 k0 p mo0s = gibbsSampler k0 mo0s [] - where (Harmonium m n) = manifold p - gibbsSampler 0 mos acc = return . reverse $ zip (potentialMapping <$> (p >$> mos)) mos:acc - gibbsSampler k mos acc = do - ls <- mapM generate $ p >$> mos - let mls = sufficientStatistic m <$> ls - os' <- mapM generate $ harmoniumTranspose p >$> mls - let mos' = sufficientStatistic n <$> os' - gibbsSampler (k-1) mos' (zip mls mos:acc) - -modulateHarmoniumBelief :: (Manifold m, Manifold n) - => Mixture :#: m - -> NaturalFunction :#: Harmonium m n - -> NaturalFunction :#: Harmonium m n --- | Adds the projection of the given belief to the biases over the state. -modulateHarmoniumBelief z trns = - let (lb,mtx,ob) = splitHarmonium trns - in joinHarmonium lb mtx $ ob <+> matrixTranspose mtx >.> z - - ---- Transducers --- - -normalBias :: (Standard :#: Normal) -> Double -normalBias sp = - let [mu,vr] = listCoordinates sp - in - mu^2/(2*vr) - -buildNormalTransducer - :: [Standard :#: Normal] -> NaturalFunction :#: Harmonium (Replicated Poisson) Normal --- | Builds a Transducer (i.e. Population Code) which is a 'Harmonium' with --- a 'Replicated' 'Poisson' latent 'Manifold'. Here the observable 'Normal' --- is 'Normal'. -buildNormalTransducer sps = - let nps = chart Natural . transition <$> sps - rp = Replicated Poisson $ length nps - lb = fromList rp $ normalBias <$> sps - ob = fromList Normal $ replicate 2 0 - tns = fromCoordinates (Tensor rp Normal) . C.concat $ coordinates <$> nps - in joinHarmonium lb tns ob - -buildReplicatedNormalTransducer - :: [Standard :#: Replicated Normal] -> NaturalFunction :#: Harmonium (Replicated Poisson) (Replicated Normal) --- | Builds a Transducer (i.e. Population Code) which is a 'Harmonium' with --- a 'Replicated' 'Poisson' latent 'Manifold'. Here the observable 'Normal' --- is 'Replicated' 'Normal'. -buildReplicatedNormalTransducer sps = - let nps = chart Natural . transition <$> sps - m = manifold $ head sps - rp = Replicated Poisson $ length nps - lb = fromList rp $ sum . mapReplicated normalBias <$> sps - ob = fromList m $ replicate (dimension m) 0 - tns = fromCoordinates (Tensor rp m) . C.concat $ coordinates <$> nps - in joinHarmonium lb tns ob - -modulateTransducerGain :: Manifold n - => Double - -> NaturalFunction :#: Harmonium (Replicated Poisson) n - -> NaturalFunction :#: Harmonium (Replicated Poisson) n --- | Multiplies the current gain of the transducer by the given value. --- Transducers are intially constructed with a gain of 1, and so initially --- this will simply set the gain. -modulateTransducerGain gn trns = - let (lb,mtx,ob) = splitHarmonium trns - lb' = alterCoordinates (+ log gn) lb - in joinHarmonium lb' mtx ob - - ---- Instances --- - - --- Harmoniums -- - -instance (Manifold m, Manifold n) => Manifold (Harmonium m n) where - dimension (Harmonium m n) = dimension m * dimension n + dimension m + dimension n - -instance (Manifold m, Manifold n) => Map (Harmonium m n) where - type Domain (Harmonium m n) = n - domain (Harmonium _ n) = n - type Codomain (Harmonium m n) = m - codomain (Harmonium m _) = m - -instance (Manifold m, Manifold n) => Apply c d (Harmonium m n) where - (>.>) p x = - let (lb,mtxp,_) = splitHarmonium p - in lb <+> (mtxp >.> x) - (>$>) p xs = - let (lb,mtxp,_) = splitHarmonium p - in (lb <+>) <$> (mtxp >$> xs) diff --git a/Goal/Probability/Graphical/NeuralNetwork.hs b/Goal/Probability/Graphical/NeuralNetwork.hs deleted file mode 100644 --- a/Goal/Probability/Graphical/NeuralNetwork.hs +++ /dev/null @@ -1,239 +0,0 @@ --- | Multilayer perceptrons and backpropagation. -module Goal.Probability.Graphical.NeuralNetwork where - - ---- Imports --- - - --- Goal -- - -import Goal.Geometry -import Goal.Probability.ExponentialFamily -import Goal.Probability.Graphical - -import qualified Data.Vector.Storable as C - - ---- Neural Networks --- - - --- | A mutlilayer perceptron with three layers. -data NeuralNetwork m n o = NeuralNetwork m n o deriving (Eq, Read, Show) - - ---- Functions --- - -splitNeuralNetwork - :: (Manifold m, Manifold n, Manifold o) - => Function Mixture Mixture :#: NeuralNetwork m n o - -> (Natural :#: m, NaturalFunction :#: Tensor m n, Natural :#: n, NaturalFunction :#: Tensor n o) --- | Splits the 'NeuralNetwork' into its component affine transformations. -splitNeuralNetwork nnp = - let (NeuralNetwork m n o) = manifold nnp - tns1 = Tensor m n - tns2 = Tensor n o - css = coordinates nnp - (mcs,css') = C.splitAt (dimension m) css - (mtx1cs,css'') = C.splitAt (dimension tns1) css' - (ncs,mtx2cs) = C.splitAt (dimension n) css'' - mp = fromCoordinates m mcs - mtx1 = fromCoordinates tns1 mtx1cs - np = fromCoordinates n ncs - mtx2 = fromCoordinates tns2 mtx2cs - in (mp,mtx1,np,mtx2) - -joinNeuralNetwork - :: (Manifold m, Manifold n, Manifold o) - => Natural :#: m - -> NaturalFunction :#: Tensor m n - -> Natural :#: n - -> NaturalFunction :#: Tensor n o - -> Function Mixture Mixture :#: NeuralNetwork m n o --- | Construct a 'NeuralNetwork' from component affine transformations. -joinNeuralNetwork mp mtx1 np mtx2 = - let (Tensor m n) = manifold mtx1 - (Tensor _ o) = manifold mtx2 - in fromCoordinates (NeuralNetwork m n o) $ coordinates mp C.++ coordinates mtx1 C.++ coordinates np C.++ coordinates mtx2 - -feedForward - :: (ExponentialFamily m, ExponentialFamily n, Manifold o) - => Function Mixture Mixture :#: NeuralNetwork m n o - -> [Mixture :#: o] - -> ([Natural :#: n], [Mixture :#: n], [Natural :#: m], [Mixture :#: m]) --- | Feeds an input forward through the network, and returns every step of --- the computation. -feedForward nnp xps = - let (mp,mtx1,np,mtx2) = splitNeuralNetwork nnp - nyps = map (<+> np) $ mtx2 >$> xps - yps = potentialMapping <$> nyps - nzps = map (<+> mp) $ mtx1 >$> yps - zps = potentialMapping <$> nzps - in (nyps,yps,nzps,zps) - -feedBackward - :: (Legendre Natural m, Legendre Natural n, Riemannian Natural m, Riemannian Natural n, Manifold o) - => Function Mixture Mixture :#: NeuralNetwork m n o - -> [Mixture :#: o] - -> [Natural :#: n] - -> [Mixture :#: n] - -> [Natural :#: m] - -> [Natural :#: m] - -> Differentials :#: Tangent (Function Mixture Mixture) (NeuralNetwork m n o) --- | Given the results of a feed forward application, back propagates a --- given error (last input) through the network. -feedBackward nnp xps nyps yps nzps errs1 = - let (_,mtx1,_,_) = splitNeuralNetwork nnp - dmps = zipWith legendreFlat nzps errs1 - dmtx1s = [ dmp >.< yp | (dmp,yp) <- zip dmps yps ] - errs2 = matrixTranspose mtx1 >$> dmps - dnps = zipWith legendreFlat nyps errs2 - dmtx2s = [ dnp >.< xp | (dnp,xp) <- zip dnps xps ] - in fromCoordinates (Tangent nnp) $ coordinates (meanPoint dmps) C.++ coordinates (meanPoint dmtx1s) - C.++ coordinates (meanPoint dnps) C.++ coordinates (meanPoint dmtx2s) - -meanSquaredBackpropagation - :: (Riemannian Natural m, Riemannian Natural n, ExponentialFamily m, ExponentialFamily n, Manifold o) - => Function Mixture Mixture :#: NeuralNetwork m n o - -> [Mixture :#: o] - -> [Mixture :#: m] - -> Differentials :#: Tangent (Function Mixture Mixture) (NeuralNetwork m n o) --- | Backpropagation algorithm with the mean squared error function. -meanSquaredBackpropagation nnp xps tps = - let (nyps,yps,nzps,zps) = feedForward nnp xps - errs1 = [ alterChart Natural $ zp <-> tp | (tp,zp) <- zip tps zps ] - in feedBackward nnp xps nyps yps nzps errs1 - - ---- Instances --- - - -instance (Manifold m, Manifold n, Manifold o) => Manifold (NeuralNetwork m n o) where - dimension (NeuralNetwork m n o) = dimension m + dimension m * dimension n + dimension n + dimension n * dimension o - -instance (ExponentialFamily m, ExponentialFamily n, Manifold o) => Map (NeuralNetwork m n o) where - type Domain (NeuralNetwork m n o) = o - domain (NeuralNetwork _ _ o) = o - type Codomain (NeuralNetwork m n o) = m - codomain (NeuralNetwork m _ _) = m - -instance (ExponentialFamily m, ExponentialFamily n, Manifold o) => Apply Mixture Mixture (NeuralNetwork m n o) where - (>$>) nnp xps = - let (_,_,_,zps) = feedForward nnp xps - in zps - - ---- Backprop --- - - -{- ---backpropagation :: NeuralNetwork (m ': ms) -> (Mixture :#: m -> Mixture :#: m) -> Differential :#: -backpropagate :: NeuralNetwork (m ': ms) -> Mixture :#: m -> Differential :#: NeuralNetwork (m ': ms) -backpropagate nnp dp = - - - ---- Internal --- - - -popManifold :: NeuralNetwork (m ': ms) -> (m, NeuralNetwork ms) -popManifold (Layer m ms) = (m,ms) - -popNeuralNetwork - :: (Manifold m, Manifold n, Manifold (NeuralNetwork (n ': ms))) - => Function Mixture Mixture :#: NeuralNetwork (m ': n ': ms) - -> (Natural :#: m, NaturalFunction :#: Tensor m n, Function Mixture Mixture :#: NeuralNetwork (n ': ms)) -popNeuralNetwork nnp = - let (m,nn') = popManifold $ manifold nnp - (n,_) = popManifold nn' - tns = Tensor m n - css = coordinates nnp - (mcs,css') = C.splitAt (dimension m) css - (mtxcs,nncs') = C.splitAt (dimension tns) css' - mp = fromCoordinates m mcs - mtx = fromCoordinates tns mtxcs - nnp' = fromCoordinates nn' nncs' - in (mp,mtx,nnp') - -feedForward - :: Function Mixture Mixture :#: NeuralNetwork ms - -> [Mixture :#: Domain (NeuralNetwork ms)] - -> [Mixture :#: Responses ms] -feedForward nnp0 xps0 = - recurse nnp0 xps0 [ chart Mixture . fromCoordinates (Responses $ Layer (manifold xp) Nub) | xp <- xps ] - where recurse nnp xps rss = - let (b,mtx,nnp') = popNeuralNetwork nnp - yps = nnp' >$> xps - in map (potentialMapping . (<+> b)) $ mtx >$> ys - - -feedBackward - :: [Mixture :#: Codomain (NeuralNetwork ms)] - -> [Mixture :#: Responses ms] - -> Differential :#: Tangent (Function Mixture Mixture) (NeuralNetwork ms) -feedBackward = undefined - ---- Instances --- - - --- Responses -- - -instance Eq (Responses '[]) where - (==) _ _ = True - -instance (Eq m, Eq (NeuralNetwork ms)) => Eq (Responses (m ': ms)) where - (==) (Responses (Layer m ms)) (Responses (Layer m' ms')) - | m == m' = ms == ms' - | otherwise = False - -instance Manifold (Responses '[]) where - dimension _ = 0 - - -instance (Manifold m, Manifold (NeuralNetwork ms)) => Manifold (Responses (m ': ms)) where - dimension (Responses (Layer m ms)) = dimension m + dimension ms - - --- NeuralNetwork -- - -instance Eq (NeuralNetwork '[]) where - (==) _ _ = True - -instance (Eq m, Eq (NeuralNetwork ms)) => Eq (NeuralNetwork (m ': ms)) where - (==) (Layer m ms) (Layer m' ms') - | m == m' = ms == ms' - | otherwise = False - -instance Manifold (NeuralNetwork '[]) where - dimension _ = 0 - -instance Manifold m => Manifold (NeuralNetwork '[m]) where - dimension _ = 0 - -instance (Manifold m, Manifold n, Manifold (NeuralNetwork (n ': ms))) => Manifold (NeuralNetwork (m ': n ': ms)) where - dimension (Layer m (Layer n ms)) = dimension m + dimension m * dimension n + dimension (Layer n ms) - -instance Manifold m => Map (NeuralNetwork '[m]) where - type Domain (NeuralNetwork '[m]) = m - domain (Layer m _) = m - type Codomain (NeuralNetwork '[m]) = m - codomain (Layer m _) = m - -instance (ExponentialFamily m, Manifold n) => Apply Mixture Mixture (NeuralNetwork '[m,n]) where - (>$>) p xs = - let (b,mtx,_) = popNeuralNetwork p - in map (potentialMapping . (<+> b)) $ mtx >$> xs - -instance (ExponentialFamily m, Manifold n, Map (NeuralNetwork (n ': ms))) - => Map (NeuralNetwork (m ': n ': ms)) where - type Domain (NeuralNetwork (m ': n ': ms)) = Domain (NeuralNetwork (n ': ms)) - domain (Layer _ nn) = domain nn - type Codomain (NeuralNetwork (m ': n ': ms)) = m - codomain (Layer m _) = m - -instance (ExponentialFamily m, Manifold n, Apply Mixture Mixture (NeuralNetwork (n ': o ': ms))) - => Apply Mixture Mixture (NeuralNetwork (m ': n ': o ': ms)) where - (>$>) p xs = - let (b,mtx,p') = popNeuralNetwork p - ys = p' >$> xs - in map (potentialMapping . (<+> b)) $ mtx >$> ys - -} diff --git a/Goal/Probability/Statistical.hs b/Goal/Probability/Statistical.hs --- a/Goal/Probability/Statistical.hs +++ b/Goal/Probability/Statistical.hs @@ -1,17 +1,28 @@ -module Goal.Probability.Statistical ( - -- * Stastical Manifolds - Statistical (sampleSpace) +{-# LANGUAGE UndecidableInstances #-} + +-- | Core types, classes, and functions for working with manifolds of +-- probability distributions. +module Goal.Probability.Statistical + ( -- * Random + Random (Random) + , Statistical (SamplePoint) , Sample - , samples - , SampleSpace - -- ** Standard Chart - , Standard (Standard) - , standardGenerate - -- ** Distributions - , Generative (generate) - , AbsolutelyContinuous (density) + , realize + -- * Initializiation + , initialize + , uniformInitialize + , uniformInitialize' + -- * Properties of Distributions + , Generative (sample,samplePoint) + , AbsolutelyContinuous (densities,logDensities) + , density + , logDensity + , Discrete (Cardinality,sampleSpace) + , pointSampleSpace , expectation + -- ** Maximum Likelihood Estimation , MaximumLikelihood (mle) + , LogLikelihood (logLikelihood,logLikelihoodDifferential) ) where @@ -20,112 +31,195 @@ -- Package -- +import Goal.Core import Goal.Geometry --- Unqualified -- +import qualified Goal.Core.Vector.Boxed as B +import qualified Goal.Core.Vector.Storable as S +import qualified Goal.Core.Vector.Generic as G -import System.Random.MWC.Monad +-- Qualified -- +import qualified Data.List as L +import qualified System.Random.MWC as R ---- Test Bed --- +import Foreign.Storable --- Probability Theory --- --- | A 'Statistical' 'Manifold' is a 'Manifold' of probability distributions, --- which all have in common a particular 'SampleSpace'. -class (Set (SampleSpace m), Manifold m) => Statistical m where - type SampleSpace m :: * - sampleSpace :: m -> SampleSpace m +-- | A 'Manifold' is 'Statistical' if it is a set of probability distributions +-- over a particular sample space, where the sample space is a set of the +-- specified 'SamplePoint's. +class Manifold x => Statistical x where + type SamplePoint x :: Type --- | A 'Sample' is an 'Element' of the 'SampleSpace'. -type Sample m = Element (SampleSpace m) +-- | A 'Sample' is a list of 'SamplePoint's. +type Sample x = [SamplePoint x] -samples :: (Discrete (SampleSpace m), Statistical m) => m -> [Sample m] --- | The list of 'Sample's. -samples = elements . sampleSpace +-- | A random variable. +newtype Random a = Random (forall s. R.Gen s -> ST s a) --- | A distribution is 'Generative' if we can 'generate' samples from it. Generation is --- powered by MWC Monad. -class Statistical m => Generative c m where - generate :: c :#: m -> RandST r (Sample m) +-- | Turn a random variable into an IO action. +realize :: Random a -> IO a +realize (Random rv) = R.withSystemRandomST rv --- | If a distribution is 'AbsolutelyContinuous' with respect to a reference --- measure on its 'SampleSpace', then we may define the 'density' of a --- probability distribution as the Radon-Nikodym derivative of the probability --- measure with respect to the base measure. -class Statistical m => AbsolutelyContinuous c m where - density :: c :#: m -> Sample m -> Double +-- | Probability distributions for which the sample space is countable. This +-- affords brute force computation of expectations. +class KnownNat (Cardinality x) => Discrete x where + type Cardinality x :: Nat + sampleSpace :: Proxy x -> Sample x --- | 'expectation' computes the brute force expected value of a 'Discrete' set given an appropriate 'density'. -expectation :: (AbsolutelyContinuous c m, Discrete (SampleSpace m)) => c :#: m -> (Sample m -> Double) -> Double +-- | Convenience function for getting the sample space of a 'Discrete' +-- probability distribution. +pointSampleSpace :: forall c x . Discrete x => c # x -> Sample x +pointSampleSpace _ = sampleSpace (Proxy :: Proxy x) + +-- | A distribution is 'Generative' if we can 'sample' from it. Generation is +-- powered by @mwc-random@. +class Statistical x => Generative c x where + samplePoint :: Point c x -> Random (SamplePoint x) + samplePoint = fmap head . sample 1 + sample :: Int -> Point c x -> Random (Sample x) + sample n = replicateM n . samplePoint + + +-- | The distributions \(P \in \mathcal M\) in a 'Statistical' 'Manifold' +-- \(\mathcal M\) are 'AbsolutelyContinuous' if there is a reference measure +-- \(\mu\) and a function \(p\) such that +-- \(P(A) = \int_A p d\mu\). We refer to \(p(x)\) as the 'density' of the +-- probability distribution. +class Statistical x => AbsolutelyContinuous c x where + logDensities :: Point c x -> Sample x -> [Double] + logDensities p = map log . densities p + + densities :: Point c x -> Sample x -> [Double] + densities p = map exp . logDensities p + +logDensity :: AbsolutelyContinuous c x => Point c x -> SamplePoint x -> Double +logDensity p = head . logDensities p . (:[]) + +density :: AbsolutelyContinuous c x => Point c x -> SamplePoint x -> Double +density p = exp . logDensity p + +-- | 'expectation' computes the brute force expected value of a 'Finite' set +-- given an appropriate 'density'. +expectation + :: forall c x . (AbsolutelyContinuous c x, Discrete x) + => Point c x + -> (SamplePoint x -> Double) + -> Double expectation p f = - let xs = elements . sampleSpace $ manifold p - in sum $ zipWith (*) (f <$> xs) (density p <$> xs) + let xs = sampleSpace (Proxy :: Proxy x) + in sum $ zipWith (*) (f <$> xs) (densities p xs) +-- Maximum Likelihood Estimation -- | 'mle' computes the 'MaximumLikelihood' estimator. -class Statistical m => MaximumLikelihood c m where - mle :: m -> [Sample m] -> c :#: m +class Statistical x => MaximumLikelihood c x where + mle :: Sample x -> c # x --- Standard Chart -- +-- | Average log-likelihood and the differential for gradient ascent. +class Manifold x => LogLikelihood c x s where + logLikelihood :: [s] -> c # x -> Double + --logLikelihood xs p = average $ log <$> densities p xs + logLikelihoodDifferential :: [s] -> c # x -> c #* x --- | A parameterization which represents the standard or typical parameterization of --- the given manifold, e.g. the 'Poisson' rate or 'Normal' mean and standard deviation. -data Standard = Standard deriving (Eq, Read, Show) -standardGenerate :: (Manifold m, Generative Standard m, Transition c Standard m) => c :#: m -> RandST r (Sample m) -standardGenerate = generate . chart Standard . transition +--- Construction --- + +-- | Generates a random point on the target 'Manifold' by generating random +-- samples from the given distribution. +initialize + :: (Manifold x, Generative d y, SamplePoint y ~ Double) + => d # y + -> Random (c # x) +initialize q = Point <$> S.replicateM (samplePoint q) + +-- | Generates an initial point on the target 'Manifold' by generating uniform +-- samples from the given vector of bounds. +uniformInitialize' :: Manifold x => B.Vector (Dimension x) (Double,Double) -> Random (Point c x) +uniformInitialize' bnds = + Random $ \gn -> Point . G.convert <$> mapM (`R.uniformR` gn) bnds + +-- | Generates an initial point on the target 'Manifold' by generating uniform +-- samples from the given vector of bounds. +uniformInitialize :: Manifold x => (Double,Double) -> Random (Point c x) +uniformInitialize bnds = + Random $ \gn -> Point <$> S.replicateM (R.uniformR bnds gn) + + --- Instances --- --- DirectSums -- +-- Random -- -instance (Statistical m, Statistical n) => Statistical (m,n) where - type SampleSpace (m,n) = (SampleSpace m, SampleSpace n) - sampleSpace (m,n) = (sampleSpace m,sampleSpace n) +instance Functor Random where + fmap f (Random rx) = + Random $ fmap f . rx -instance (Generative c m, Generative c n) => Generative c (m,n) where - generate cmn = do - let (cm,cn) = splitPair' cmn - mx <- generate cm - nx <- generate cn - return (mx, nx) +instance Applicative Random where + pure x = Random $ \_ -> return x + (<*>) = ap -instance (AbsolutelyContinuous Standard m, AbsolutelyContinuous Standard n) => AbsolutelyContinuous Standard (m,n) where - density cmn (mx,nx) = - let (cm,cn) = splitPair' cmn - in density cm mx * density cn nx +instance Monad Random where + (>>=) (Random rx) rf = + Random $ \gn -> do + a <- rx gn + let (Random rv) = rf a + rv gn + -- Replicated -- -instance Statistical m => Statistical (Replicated m) where - type SampleSpace (Replicated m) = Replicated (SampleSpace m) - sampleSpace (Replicated m n) = Replicated (sampleSpace m) n +instance (Statistical x, KnownNat k, Storable (SamplePoint x)) + => Statistical (Replicated k x) where + type SamplePoint (Replicated k x) = S.Vector k (SamplePoint x) -instance (Statistical m, Generative c m) => Generative c (Replicated m) where - generate = sequence . mapReplicated generate +instance (KnownNat k, Generative c x, Storable (SamplePoint x)) + => Generative c (Replicated k x) where + samplePoint = S.mapM samplePoint . splitReplicated -instance (Statistical m, AbsolutelyContinuous Standard m) => AbsolutelyContinuous Standard (Replicated m) where - density ds xs = product $ zipWith ($) (mapReplicated density ds) xs +instance (KnownNat k, Storable (SamplePoint x), AbsolutelyContinuous c x) + => AbsolutelyContinuous c (Replicated k x) where + densities cx sxss = + let sxss' = L.transpose $ S.toList <$> sxss + cxs = S.toList $ splitReplicated cx + dnss = zipWith densities cxs sxss' + in product <$> L.transpose dnss -instance (Statistical m, Transition Standard c m) => Transition Standard c (Replicated m) where - transition = joinReplicated . mapReplicated transition +instance (KnownNat k, LogLikelihood c x s, Storable s) + => LogLikelihood c (Replicated k x) (S.Vector k s) where + logLikelihood cxs ps = S.sum . S.imap subLogLikelihood $ splitReplicated ps + where subLogLikelihood fn = logLikelihood (flip S.index fn <$> cxs) + logLikelihoodDifferential cxs ps = + joinReplicated . S.imap subLogLikelihoodDifferential $ splitReplicated ps + where subLogLikelihoodDifferential fn = + logLikelihoodDifferential (flip S.index fn <$> cxs) -instance (Statistical m, Transition c Standard m) => Transition c Standard (Replicated m) where - transition = joinReplicated . mapReplicated transition +-- Pair -- ---- Graveyard --- +instance (Statistical x) => Statistical [x] where + type SamplePoint [x] = [SamplePoint x] -{- -manifoldExpectation :: (Manifold n, AbsolutelyContinuous c m, Discrete (SampleSpace m)) - => c :#: m -> (Sample m -> d :#: n) -> d :#: n -manifoldExpectation p f = - let xs = elements . sampleSpace $ manifold p - in foldl1' (<+>) $ zipWith (.>) (density p <$> xs) (f <$> xs) +instance (Statistical x, Statistical y) + => Statistical (x,y) where + type SamplePoint (x,y) = (SamplePoint x, SamplePoint y) --} +instance (Generative c x, Generative c y) => Generative c (x,y) where + samplePoint pmn = do + let (pm,pn) = split pmn + xm <- samplePoint pm + xn <- samplePoint pn + return (xm,xn) + +instance (AbsolutelyContinuous c x, AbsolutelyContinuous c y) + => AbsolutelyContinuous c (x,y) where + densities cxy sxys = + let (cx,cy) = split cxy + (sxs,sys) = unzip sxys + in zipWith (*) (densities cx sxs) $ densities cy sys diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,72 @@ +This library provides tools for implementing and applying statistical and +machine learning algorithms. The core concept of goal-probability is that of a +statistical manifold, i.e. manifold of probability distributions, with a focus +on exponential family distributions. What follows is brief introduction to how +we define and work with statistical manifolds in Goal. + +The core definition of this library is that of a `Statistical` `Manifold`. +```haskell +class Manifold x => Statistical x where + type SamplePoint x :: Type +``` +A `Statistical` `Manifold` is a `Manifold` of probability distributions, such +that each point on the manifold is a probability distribution over the specified +space of `SamplePoint`s. We may evaluate the probability mass/density of a `SamplePoint` under a given distribution as long as the distribution is `AbsolutelyContinous`. +```haskell +class Statistical x => AbsolutelyContinuous c x where + density :: Point c x -> SamplePoint x -> Double + densities :: Point c x -> Sample x -> [Double] +``` +Similarly, we may generate a `Sample` from a given distribution as long as it is `Generative`. +```haskell +type Sample x = [SamplePoint x] + +class Statistical x => Generative c x where + samplePoint :: Point c x -> Random r (SamplePoint x) + sample :: Int -> Point c x -> Random r (Sample x) +``` +In both these cases, class methods are defined both both single and bulk +evaluation, to potentially take advantage of bulk linear algebra operations. + +Let us now look at some example distributions that we may define; for the sake +of brevity, I will not introduce every bit of necessary code. In +Goal we create a normal distribution by writing +```haskell +nrm :: Source # Normal +nrm = fromTuple (0,1) +``` +where 0 is the mean and 1 is the variance. For each `Statistical` `Manifold`, +the `Source` coordinate system represents some standard parameterization, e.g. +as one typically finds on Wikipedia. Similarly, we can create a binomial +distribution with +```haskell +bnm :: Source # Binomial 5 +bnm = Point $ S.singleton 0.5 +``` +which is a binomial distribution over 5 fair coin tosses. + +Exponential families are also a core part of this library. An `ExponentiaFamily` +is a kind of `Statistical` `Manifold` defined in terms of a +`sufficientStatistic` and a `baseMeasure`. +```haskell +class Statistical x => ExponentialFamily x where + sufficientStatistic :: SamplePoint x -> Mean # x + baseMeasure :: Proxy x -> SamplePoint x -> Double +``` + +Exponential families may always be parameterized in terms of the so-called +`Natural` and `Mean` parameters. Mean coordinates are equal to the average value +of the `sufficientStatistic` under the given distribution. The `Natural` +coordinates are arguably less intuitive, but they are critical for implementing +evaluating exponential family distributions numerically. For example, the +unnormalized density function of an `ExponentialFamily` distribution is +given by +```haskell +unnormalizedDensity :: forall x . ExponentialFamily x => Natural # x -> SamplePoint x -> Double +unnormalizedDensity p x = + exp (p <.> sufficientStatistic x) * baseMeasure (Proxy @ x) x +``` + +For in-depth tutorials visit my +[blog](https://sacha-sokoloski.gitlab.io/website/pages/blog.html). + diff --git a/benchmarks/backpropagation.hs b/benchmarks/backpropagation.hs new file mode 100644 --- /dev/null +++ b/benchmarks/backpropagation.hs @@ -0,0 +1,79 @@ +{-# LANGUAGE TypeOperators,TypeFamilies,FlexibleContexts,DataKinds #-} + +--- Imports --- + + +-- Goal -- + +import Goal.Core +import Goal.Geometry +import Goal.Probability + +import qualified Goal.Core.Vector.Storable as S + +-- Qualified -- + +import qualified Criterion.Main as C + + +--- Globals --- + + +-- Data -- + +f :: Double -> Double +f x = exp . sin $ 2 * x + +mnx,mxx :: Double +mnx = -3 +mxx = 3 + +xs :: [Double] +xs = range mnx mxx 200 + +fp :: Source # Normal +fp = Point $ S.doubleton 0 0.1 + +-- Neural Network -- + +cp :: Source # Normal +cp = Point $ S.doubleton 0 0.0001 + +type NeuralNetwork' = + NeuralNetwork '[ '(Tensor, R 1000 Bernoulli), '(Tensor, R 1000 Bernoulli)] + Tensor NormalMean NormalMean + + + +-- Training -- + +nepchs :: Int +nepchs = 1 + +eps :: Double +eps = 0.0001 + +-- Layout -- + +main :: IO () +main = do + + ys <- realize $ mapM (noisyFunction fp f) xs + + mlp0 <- realize $ initialize cp + + let xys = zip ys xs + + let cost :: Natural # NeuralNetwork' -> Double + cost = conditionalLogLikelihood xys + + let backprop :: Natural # NeuralNetwork' -> Natural #* NeuralNetwork' + backprop = conditionalLogLikelihoodDifferential xys + + admmlps0 mlp = take nepchs $ vanillaGradientSequence backprop eps defaultAdamPursuit mlp + + let mlp = last $!! admmlps0 mlp0 + + C.defaultMain + [ C.bench "application" $ C.nf cost mlp + , C.bench "backpropagation" $ C.nf backprop mlp ] diff --git a/benchmarks/regression.hs b/benchmarks/regression.hs new file mode 100644 --- /dev/null +++ b/benchmarks/regression.hs @@ -0,0 +1,112 @@ +{-# LANGUAGE TypeOperators,TypeFamilies,FlexibleContexts,DataKinds #-} + + +--- Imports --- + + +-- Goal -- + +import Goal.Core +import Goal.Geometry +import Goal.Probability + +import qualified Goal.Core.Vector.Storable as S + +-- Qualified -- + +import qualified Criterion.Main as C + + +--- Globals --- + + +-- Data -- + +f :: Double -> Double +f x = exp . sin $ 2 * x + +mnx,mxx :: Double +mnx = -3 +mxx = 3 + +xs :: [Double] +xs = concat . replicate 5 $ range mnx mxx 8 + +fp :: Source # Normal +fp = Point $ S.doubleton 0 0.1 + +-- Neural Network -- + +cp :: Source # Normal +cp = Point $ S.doubleton 0 0.1 + +type NeuralNetwork' = + NeuralNetwork '[ '(Tensor, R 50 Bernoulli)] + Tensor NormalMean NormalMean + +-- Training -- + +nepchs :: Int +nepchs = 1000 + +eps :: Double +eps = 0.01 + +-- Momentum +mxmu :: Double +mxmu = 0.999 + + +--- Main --- + + +main :: IO () +main = do + + ys <- realize $ mapM (noisyFunction fp f) xs + + mlp0 <- realize $ initialize cp + + let xys = zip ys xs + + let cost :: Natural # NeuralNetwork' -> Double + cost = conditionalLogLikelihood xys + + let backprop :: Natural # NeuralNetwork' -> Natural #* NeuralNetwork' + backprop = conditionalLogLikelihoodDifferential xys + + let sortedBackprop :: Natural # NeuralNetwork' -> Natural #* NeuralNetwork' + sortedBackprop = mapConditionalLogLikelihoodDifferential $ conditionalDataMap xys + + sgdmlps0 mlp = take nepchs $ mlp0 : vanillaGradientSequence backprop eps Classic mlp + mtmmlps0 mlp = take nepchs + $ mlp0 : vanillaGradientSequence backprop eps (defaultMomentumPursuit mxmu) mlp + admmlps0 mlp = take nepchs + $ mlp0 : vanillaGradientSequence backprop eps defaultAdamPursuit mlp + sadmmlps0 mlp = take nepchs + $ mlp0 : vanillaGradientSequence sortedBackprop eps defaultAdamPursuit mlp + + C.defaultMain + [ C.bench "sgd" $ C.nf sgdmlps0 mlp0 + , C.bench "momentum" $ C.nf mtmmlps0 mlp0 + , C.bench "adam" $ C.nf admmlps0 mlp0 + , C.bench "sorted-adam" $ C.nf sadmmlps0 mlp0 ] + + let sgdmlps = sgdmlps0 mlp0 + mtmmlps = mtmmlps0 mlp0 + admmlps = admmlps0 mlp0 + sadmmlps = sadmmlps0 mlp0 + + let sgdcst = cost $ last sgdmlps + mtmcst = cost $ last mtmmlps + admcst = cost $ last admmlps + sadmcst = cost $ last sadmmlps + + putStrLn "SGD LL:" + print sgdcst + putStrLn "Momentum LL:" + print mtmcst + putStrLn "Adam LL:" + print admcst + putStrLn "Sorted Adam LL:" + print sadmcst diff --git a/goal-probability.cabal b/goal-probability.cabal --- a/goal-probability.cabal +++ b/goal-probability.cabal @@ -1,134 +1,81 @@ +cabal-version: 3.0 +version: 0.20 name: goal-probability -version: 0.1 -synopsis: Manifolds of probability distributions -description: Provides probability distributions, exponential families, as well - as things based on exponential families such as multilayer perceptrons and - harmoniums (e.g. restricted Boltzmann machines). -license: BSD3 +synopsis: Optimization on manifolds of probability distributions with Goal +description: goal-probability provides tools for implementing and applying basic statistical models. The core concept of goal-probability are statistical manifolds, i.e. manifold of probability distributions, with a focus on exponential family distributions. +license: BSD-3-Clause license-file: LICENSE +extra-source-files: README.md author: Sacha Sokoloski -maintainer: sokolo@mis.mpg.de +maintainer: sacha.sokoloski@mailbox.org +homepage: https://gitlab.com/sacha-sokoloski/goal category: Math build-type: Simple -cabal-version: >=1.10 library exposed-modules: - Goal.Probability, - Goal.Probability.Distributions, - Goal.Probability.ExponentialFamily, - Goal.Probability.Statistical, - Goal.Probability.Graphical, - Goal.Probability.Graphical.Harmonium, - Goal.Probability.Graphical.NeuralNetwork - default-extensions: TypeOperators, TypeFamilies, FlexibleInstances, - FlexibleContexts, MultiParamTypeClasses - build-depends: - base==4.*, - mwc-random==0.13.*, - mwc-random-monad==0.7.*, - math-functions==0.1.5.*, - vector==0.11.*, - hmatrix==0.17.*, - statistics==0.13.*, - goal-core==0.1, - goal-geometry==0.1 - default-language: Haskell2010 - ghc-options: -O2 -Wall -fno-warn-type-defaults -fno-warn-missing-signatures - -executable cross-entropy-descent - main-is: cross-entropy-descent.hs - hs-source-dirs: scripts - ghc-options: -Wall -O2 -threaded -rtsopts -fno-warn-type-defaults - -fno-warn-missing-signatures -fno-warn-unused-do-bind - build-depends: - base==4.*, - goal-core==0.1, - goal-geometry==0.1, - goal-probability==0.1 - default-language: Haskell2010 - -executable poisson-binomial - main-is: poisson-binomial.hs - hs-source-dirs: scripts - ghc-options: -Wall -O2 -threaded -rtsopts -fno-warn-type-defaults - -fno-warn-missing-signatures -fno-warn-unused-do-bind - build-depends: - base==4.*, - goal-core==0.1, - goal-geometry==0.1, - goal-probability==0.1 - default-language: Haskell2010 - -executable univariate - main-is: univariate.hs - hs-source-dirs: scripts - ghc-options: -Wall -O2 -threaded -rtsopts -fno-warn-type-defaults - -fno-warn-missing-signatures -fno-warn-unused-do-bind - build-depends: - base==4.*, - goal-core==0.1, - goal-geometry==0.1, - goal-probability==0.1 - default-language: Haskell2010 - -executable multivariate - main-is: multivariate.hs - hs-source-dirs: scripts - ghc-options: -Wall -O2 -threaded -rtsopts -fno-warn-type-defaults - -fno-warn-missing-signatures -fno-warn-unused-do-bind - build-depends: - base==4.*, - goal-core==0.1, - goal-geometry==0.1, - goal-probability==0.1, - vector==0.11.* - default-language: Haskell2010 - -executable transducer - main-is: transducer.hs - hs-source-dirs: scripts - ghc-options: -Wall -O2 -threaded -rtsopts -fno-warn-type-defaults - -fno-warn-missing-signatures -fno-warn-unused-do-bind - build-depends: - base==4.*, - goal-core==0.1, - goal-geometry==0.1, - goal-probability==0.1 - default-language: Haskell2010 - -executable transducer-field - main-is: transducer-field.hs - hs-source-dirs: scripts - ghc-options: -Wall -O2 -threaded -rtsopts -fno-warn-type-defaults - -fno-warn-missing-signatures -fno-warn-unused-do-bind + Goal.Probability + Goal.Probability.Statistical + Goal.Probability.ExponentialFamily + Goal.Probability.Distributions + Goal.Probability.Distributions.CoMPoisson + Goal.Probability.Distributions.Gaussian + Goal.Probability.Conditional build-depends: - base==4.*, - goal-core==0.1, - goal-geometry==0.1, - goal-probability==0.1 + base >= 4.13 && < 4.15, + mwc-random, + hmatrix-special, + ghc-typelits-knownnat, + ghc-typelits-natnormalise, + goal-core, + parallel, + statistics, + vector, + hmatrix, + containers, + goal-geometry default-language: Haskell2010 + default-extensions: + NoStarIsType, + ScopedTypeVariables, + ExplicitNamespaces, + TypeOperators, + KindSignatures, + DataKinds, + RankNTypes, + TypeFamilies, + GeneralizedNewtypeDeriving, + StandaloneDeriving, + FlexibleContexts, + MultiParamTypeClasses, + ConstraintKinds, + FlexibleInstances + ghc-options: -Wall -O2 -executable divergence - main-is: divergence.hs - hs-source-dirs: scripts - ghc-options: -Wall -O2 -threaded -rtsopts -fno-warn-type-defaults - -fno-warn-missing-signatures -fno-warn-unused-do-bind +benchmark regression + type: exitcode-stdio-1.0 + main-is: regression.hs + hs-source-dirs: benchmarks build-depends: - base==4.*, - goal-core==0.1, - goal-geometry==0.1, - goal-probability==0.1 + base, + goal-core, + goal-geometry, + goal-probability, + bytestring, + cassava, + criterion default-language: Haskell2010 + ghc-options: -Wall -O2 -executable backpropagation +benchmark backpropagation + type: exitcode-stdio-1.0 main-is: backpropagation.hs - hs-source-dirs: scripts - ghc-options: -Wall -O2 -threaded -rtsopts -fno-warn-type-defaults - -fno-warn-missing-signatures -fno-warn-unused-do-bind + hs-source-dirs: benchmarks build-depends: - base==4.*, - goal-core==0.1, - goal-geometry==0.1, - goal-probability==0.1 + base, + goal-core, + goal-geometry, + goal-probability, + criterion default-language: Haskell2010 + ghc-options: -Wall -O2 diff --git a/scripts/backpropagation.hs b/scripts/backpropagation.hs deleted file mode 100644 --- a/scripts/backpropagation.hs +++ /dev/null @@ -1,120 +0,0 @@ -{-# LANGUAGE TypeOperators, TypeFamilies, FlexibleContexts #-} - ---- Imports --- - - --- Goal -- - -import Goal.Core -import Goal.Geometry -import Goal.Probability - - ---- Globals --- - -f x = exp . sin $ 2 * x -nsmps = 20 -mnx = -3 -mxx = 3 -xs = range mnx mxx nsmps - --- Neural Network -- - -m = Poisson -n = Replicated Bernoulli 20 -o = MeanNormal 1 - -nn = NeuralNetwork m n o - --- Training -- - -eps = 0.05 -nepchs = 10000 - --- Plot -- - -nplts = 100 -pltrng = range mnx mxx nplts - --- Layout -- - -main = do - - smps <- runWithSystemRandom $ mapM (noisyFunction (chart Standard $ fromList Normal [0,0.1]) f) xs - let xps = sufficientStatistic o <$> xs - tps = [ fromList Poisson [smp] | smp <- smps ] - - cs0 <- runWithSystemRandom . replicateM (dimension nn) . generate . chart Standard $ fromList Normal [0,0.1] - let nnp0 = fromList nn cs0 - - let gradient nnp = meanSquaredBackpropagation nnp xps tps - nnps = vanillaGradientDescent eps gradient nnp0 - nnp1 = nnps !! nepchs - - fhat x = coordinate 0 $ nnp1 >.> sufficientStatistic o x - - let lyt1 = execEC $ do - - layout_title .= "Regression" - - plot . liftEC $ do - - plot_lines_title .= "True" - plot_lines_style .= solidLine 3 (opaque black) - plot_lines_values .= [zip pltrng (f <$> pltrng)] - - plot . liftEC $ do - - plot_points_title .= "Samples" - plot_points_style .= filledCircles 4 (opaque black) - plot_points_values .= zip xs smps - - plot . liftEC $ do - - plot_lines_title .= "MLP" - plot_lines_style .= solidLine 3 (opaque red) - plot_lines_values .= [zip pltrng (fhat <$> pltrng)] - - let (mp,mtx1,np,mtx2) = splitNeuralNetwork nnp1 - let lyt2 = coordinateLogHistogram 10 "Network Weights" ["B1","I1","B2","I2"] - [coordinates mp, coordinates mtx1, coordinates np, coordinates mtx2] - - renderableToAspectWindow False 800 800 . toRenderable . weights (1,1) $ tval lyt2 ./. tval lyt1 - -{- - let hstplt = histogramPlot nb mn mx [toDouble <$> smps] . execEC $ do - plot_bars_titles .= ["Samples"] - plot_bars_item_styles .= [(solidFillStyle $ opaque blue, Nothing)] - - return . histogramLayoutLR hstplt . execEC $ do - - layoutlr_title .= (show (manifold p) ++ "; KLD: " ++ take 5 (showFFloat (Just 3) (klDivergence mle1 p) "")) - layoutlr_left_axis . laxis_title .= "Sample Count" - layoutlr_right_axis . laxis_title .= "Probability Mass" - layoutlr_x_axis . laxis_title .= "Value" - - plotRight . liftEC $ do - plot_lines_style .= dashedLine 3 [2,1] (opaque black) - plot_lines_title .= "True" - plot_lines_values .= [lineFun1 p] - - plotRight . liftEC $ do - plot_lines_style .= dashedLine 3 [10,5] (opaque red) - plot_lines_title .= "Standard MLE" - plot_lines_values .= [ lineFun1 mle1 ] - - plotRight . liftEC $ do - plot_lines_style .= dashedLine 3 [7,3] (opaque purple) - plot_lines_title .= "Exponential Family MLE" - plot_lines_values .= [ lineFun2 . chart Natural $ mle m smps ] - - lytB <- tval <$> generateLayout bnsB mnB mxB toDoubleB rngB truB - lytC <- tval <$> generateLayout bnsC mnC mxC toDoubleC rngC truC - lytP <- tval <$> generateLayout bnsP mnP mxP toDoubleP rngP truP - lytN <- tval <$> generateLayout bnsN mnN mxN toDoubleN rngN truN - - let grd1 = lytB .|. lytC - grd2 = lytP .|. lytN - - renderableToAspectWindow False 800 600 . toRenderable . weights (1,1) $ grd1 ./. grd2 - -} diff --git a/scripts/cross-entropy-descent.hs b/scripts/cross-entropy-descent.hs deleted file mode 100644 --- a/scripts/cross-entropy-descent.hs +++ /dev/null @@ -1,114 +0,0 @@ -{-# LANGUAGE TypeOperators, TypeFamilies, FlexibleContexts #-} - ---- Imports --- - - --- Goal -- - -import Goal.Core -import Goal.Geometry -import Goal.Probability - - ---- Globals --- - -nsmps = 20 - --- True Normal -- - -sp = chart Standard $ fromList Normal [1.5,2] - --- Gradient Ascent -- - -eps = 0.01 -stps = 3000 -sp0 = chart Standard $ fromList Normal [0,1] - --- Plot -- - -mnmu = 0 -mxmu = 3 -mnvr = 1 -mxvr = 4 - -axprms = LinearAxisParams (show . round) 4 4 - -m1rng = (mnmu,mxmu,600) -m2rng = (mnvr,mxvr,600) -niso = 20 -clrs = rgbaGradient (0,0,0,1) (1,0,0,1) niso - --- Functions -- - -logLikelihood p xs = sum $ log . density p <$> xs - -naturalDerivatives :: [Double] -> Natural :#: Normal -> Differentials :#: Tangent Natural Normal -naturalDerivatives xs p = fromCoordinates (Tangent p) . coordinates - $ meanPoint (sufficientStatistic Normal <$> xs) <-> potentialMapping p - -standardDerivatives :: [Double] -> Standard :#: Normal -> Differentials :#: Tangent Standard Normal -standardDerivatives xs p = - let [mu,vr] = listCoordinates p - in meanPoint [ fromList (Tangent p) [ recip vr * (xi - mu), recip (2*vr) * (recip vr * (xi - mu)^2 - 1) ] | xi <- xs ] - --- Layout -- - -main = do - - smps <- runWithSystemRandom . replicateM nsmps $ generate sp - - let mp' = chart Mixture . meanPoint $ sufficientStatistic Normal <$> smps - sp' = chart Standard $ transition mp' - - let vsps1 = take stps $ vanillaGradientAscent eps (standardDerivatives smps) sp0 - nsps1 = take stps $ gradientAscent eps (standardDerivatives smps) sp0 - - let np0 = chart Natural $ transition sp0 - vnps2 = take stps $ vanillaGradientAscent eps (naturalDerivatives smps) np0 - --nnps2 = take stps $ gradientAscent eps (naturalDerivatives smps) np0 - vsps2 = chart Standard . transition <$> vnps2 - --nsps2 = chart Standard . transition <$> nnps2 - - let rnbl = toRenderable . execEC $ do - - let f x y = logLikelihood (chart Standard $ fromList Normal [x,y]) smps - cntrs = contours m1rng m2rng niso f - - layout_x_axis . laxis_generate .= scaledAxis axprms (mnmu,mxmu) - layout_x_axis . laxis_override .= axisGridHide - layout_x_axis . laxis_title .= "μ" - layout_y_axis . laxis_generate .= scaledAxis axprms (mnvr,mxvr) - layout_y_axis . laxis_override .= axisGridHide - layout_y_axis . laxis_title .= "σ^2" - - sequence_ $ do - - ((_,cntr),clr) <- zip cntrs clrs - - return . plot . liftEC $ do - - plot_lines_style .= solidLine 3 clr - plot_lines_values .= cntr - - plot . liftEC $ do - plot_lines_style .= solidLine 3 (opaque blue) - plot_lines_values .= [toPair <$> vsps2] - - plot . liftEC $ do - plot_lines_style .= solidLine 3 (opaque green) - plot_lines_values .= [toPair <$> vsps1] - - plot . liftEC $ do - plot_lines_style .= solidLine 3 (opaque purple) - plot_lines_values .= [toPair <$> nsps1] - - plot . liftEC $ do - plot_points_style .= filledCircles 4 (opaque black) - plot_points_values .= [toPair sp] - - plot . liftEC $ do - plot_points_style .= filledCircles 4 (opaque red) - plot_points_values .= [toPair sp'] - - --renderableToAspectWindow False 800 600 . toRenderable $ lyt - void $ renderableToFile (FileOptions (500,350) PDF) "cross-entropy-descent.pdf" rnbl diff --git a/scripts/divergence.hs b/scripts/divergence.hs deleted file mode 100644 --- a/scripts/divergence.hs +++ /dev/null @@ -1,81 +0,0 @@ -{-# LANGUAGE FlexibleContexts,TypeOperators #-} - ---- Imports --- - - --- Scientific -- - -import Goal.Core -import Goal.Geometry -import Goal.Probability - ---- Program --- - - --- Globals -- - -res = 200 -niso = 10 - - --- Functions -- - -divergenceLayout :: (ExponentialFamily m, Transition c Mixture m, Transition c Natural m) - => (Double, Double) -> AlphaColour Double -> c -> m -> Layout Double Double -divergenceLayout (mn,mx) clr c m = execEC $ do - - let f x y = relativeEntropy (chart c $ fromList m [x]) (chart c $ fromList m [y]) - cntrs = contours (mn,mx,res) (mn,mx,res) niso f - x0 = (mx + mn) / 2 - y0 = x0 - str0 = "0.0" - hgh = 0.95 * mx + 0.05 * mn - lw = 0.05 * mx + 0.95 * mn - x1 = hgh - y1 = lw - str1 = showFFloat (Just 1) (f x1 y1) "" - x2 = lw - y2 = hgh - str2 = showFFloat (Just 1) (f x2 y2) "" - - plot . liftEC $ do - plot_lines_style .= solidLine 2 clr - plot_lines_values .= [[ (x,x) | x <- range mn mx 3 ]] - - sequence_ $ do - - (_,cntr) <- cntrs - - return . plot . liftEC $ do - - plot_lines_style .= solidLine 3 clr - plot_lines_values .= cntr - - plot . liftEC $ do - plot_points_values .= [(x0,y0),(x1,y1),(x2,y2)] - plot_points_style .= filledCircles 9 (opaque white) - - plot . liftEC $ do - plot_annotation_values .= [(x0,y0,str0),(x1,y1,str1),(x2,y2,str2)] - plot_annotation_style . font_weight .= FontWeightBold - - --- Main -- - -main :: IO () -main = do - - let [blyt0,blyt1,plyt0,plyt1] = - [ toRenderable $ divergenceLayout (0.02,0.98) (opaque blue) Mixture Bernoulli - , toRenderable $ divergenceLayout (-5,5) (opaque red) Natural Bernoulli - , toRenderable $ divergenceLayout (0.1,4) (opaque blue) Mixture Poisson - , toRenderable $ divergenceLayout (-2,2) (opaque red) Natural Poisson ] - - let bgrd = tval blyt0 ./. tval blyt1 - pgrd = tval plyt0 ./. tval plyt1 - - let rnbl = gridToRenderable . weights (1,1) $ bgrd .|. pgrd - --void $ renderableToFile (FileOptions (500,500) PDF) "divergence.pdf" grd - void $ renderableToAspectWindow False 1000 1000 rnbl - - diff --git a/scripts/multivariate.hs b/scripts/multivariate.hs deleted file mode 100644 --- a/scripts/multivariate.hs +++ /dev/null @@ -1,100 +0,0 @@ ---- Imports --- - - --- Goal -- - -import Goal.Core -import Goal.Geometry -import Goal.Probability - -import qualified Data.Vector.Storable as C - - ---- Globals --- - - -nsmps = 10 -tru = chart Standard $ fromList (MultivariateNormal 2) [0,0.5,1,0.5,0,1] - -rng = (-4,4,400) -niso = 10 - -axprms = LinearAxisParams (show . round) 5 5 - -vectorToPair xs = (xs C.! 0, xs C.! 1) -pairToVector (x,y) = C.fromList [x,y] - ---- Main --- - - -main :: IO () -main = do - - smps <- runWithSystemRandom . replicateM nsmps $ generate tru - - let mlenrm = chart Standard $ mle (MultivariateNormal 2) smps - --efnrm = chart Natural $ mle (MultivariateNormal 2) smps - - truf x y = density tru $ pairToVector (x,y) - mlef x y = density mlenrm $ pairToVector (x,y) - --eff x y = density efnrm $ pairToVector (x,y) - - trucntrs = contours rng rng niso truf - mlecntrs = contours rng rng niso mlef - --efcntrs = contours rng rng niso eff - - truclrs = rgbaGradient (1,0,0,0.5) (1,0,0,1) niso - mleclrs = rgbaGradient (0,0,1,0.5) (0,0,1,1) niso - --efclrs = rgbaGradient (0,1,0,0.5) (0,1,0,1) niso - bls = True : repeat False - - rnbl = toRenderable . execEC $ do - - --layout_title .= ("Multivariate Normal" ++ "; KLD: " ++ showFFloat (Just 3) (klDivergence mlenrm tru) "") - - layout_x_axis . laxis_generate .= scaledAxis axprms (-4,4) - layout_x_axis . laxis_override .= axisGridHide - layout_x_axis . laxis_title .= "x" - layout_y_axis . laxis_generate .= scaledAxis axprms (-4,4) - layout_y_axis . laxis_override .= axisGridHide - layout_y_axis . laxis_title .= "y" - - sequence_ $ do - - ((_,cntr),clr,bl) <- zip3 trucntrs truclrs bls - - return . plot . liftEC $ do - - --when bl $ plot_lines_title .= "True" - plot_lines_style .= solidLine 3 clr - plot_lines_values .= cntr - - sequence_ $ do - - ((_,cntr),clr,bl) <- zip3 mlecntrs mleclrs bls - - return . plot . liftEC $ do - - --when bl $ plot_lines_title .= "Standard MLE" - plot_lines_style .= solidLine 3 clr - plot_lines_values .= cntr - - plot . liftEC $ do - --plot_points_title .= "Samples" - plot_points_values .= map vectorToPair smps - plot_points_style .= filledCircles 4 (opaque black) - -{- - sequence $ do - - ((_,cntr),clr,bl) <- zip3 efcntrs efclrs bls - - return . plot . liftEC $ do - - when bl $ plot_lines_title .= "Exponential Family MLE" - plot_lines_style .= solidLine 3 clr - plot_lines_values .= cntr - -} - - --renderableToAspectWindow False 800 600 rnbl - void $ renderableToFile (FileOptions (250,250) PDF) "multivariate.pdf" rnbl diff --git a/scripts/poisson-binomial.hs b/scripts/poisson-binomial.hs deleted file mode 100644 --- a/scripts/poisson-binomial.hs +++ /dev/null @@ -1,51 +0,0 @@ --- A script which demonstrates how the binomial and poisson distributions --- approximate each other. - ---- Imports --- - - --- Goal -- - -import Goal.Core -import Goal.Geometry -import Goal.Probability - - ---- Script --- - - -main = renderableToAspectWindow False 800 600 . toRenderable $ poissonLayout 5 - -poissonLayout :: Double -> Layout Int Double -poissonLayout lmda = execEC $ do - - layout_title .= "Binomial Convergence to Poisson" - layout_y_axis . laxis_title .= "Probability Mass" - layout_x_axis . laxis_title .= "Count" - - let rng = [0..20] - - plot . liftEC $ do - - let pd = chart Standard $ fromList Poisson [lmda] - ppnts = zip rng $ density pd <$> rng - - plot_points_style .= filledCircles 8 (opaque red) - plot_points_title .= ("λ = " ++ show lmda) - plot_points_values .= ppnts - - let bplt n = liftEC $ do - - let p = lmda / fromIntegral n - alph = 2 * fromIntegral n / 100 - - bd = chart Standard $ fromList (Binomial n) [p] - bpnts = zip rng $ density bd <$> take (n+1) rng - - plot_points_style .= filledCircles 5 (withOpacity black alph) - plot_points_title .= ("n = " ++ show n ++ ", p = " ++ show p) - plot_points_values .= bpnts - - plot $ bplt 10 - plot $ bplt 25 - plot $ bplt 100 diff --git a/scripts/transducer-field.hs b/scripts/transducer-field.hs deleted file mode 100644 --- a/scripts/transducer-field.hs +++ /dev/null @@ -1,84 +0,0 @@ ---- Imports --- - - --- Goal -- - -import Goal.Core -import Goal.Geometry -import Goal.Probability - - ---- Program --- - - --- Globals -- - -mnx = -4 -mxx = 4 -mny = -4 -mxy = 4 -vr = 2 -sps = [ joinReplicated [fromList Normal [x,vr], fromList Normal [y,vr]] - | x <- tail $ range mnx mxx 10, y <- range mny mxy 10 ] -gn = 10 -trns = modulateTransducerGain gn $ buildReplicatedNormalTransducer sps - -x0 = -1 -y0 = 1 -xy0 = [x0,y0] - --- Functions -- - -rngx = (mnx,mxx,100) -rngy = (mny,mxy,100) -niso = 10 -clrs = rgbaGradient (0,0,1,0.6) (1,0,0,0.6) niso - -transducerRenderable rs = toRenderable . execEC $ do - - let [x',_,y',_] = listCoordinates $ conditionalObservableDistribution trns rs - posterior x y = density (conditionalObservableDistribution trns rs) [x,y] - cntrs = contours rngx rngy niso posterior - - sequence_ $ do - - ((_,cntr),clr) <- zip cntrs clrs - - return . plot . liftEC $ do - - plot_lines_style .= solidLine 3 clr - plot_lines_values .= cntr - - layout_x_axis . laxis_generate .= scaledAxis def (mnx,mxx) - layout_y_axis . laxis_generate .= scaledAxis def (mny,mxy) - - plot . liftEC $ do - plot_points_style .= filledCircles 4 (opaque black) - plot_points_values .= [(x0, y0)] - plot_points_title .= "Stimulus" - - plot . liftEC $ do - plot_points_style .= filledCircles 4 (opaque red) - plot_points_values .= [(x',y')] - plot_points_title .= "Estimate" - - plot . liftEC $ - plot_annotation_values .= [(x,y,show r) | (r,[x,_,y,_]) <- zip rs $ listCoordinates <$> sps ] - -{- - plotLeft . liftEC $ do - plot_lines_style .= solidLine 3 (opaque red) - plot_lines_values .= [let plts = posterior <$> pltrng in zip pltrng $ (*50) . (/ sum plts) <$> plts ] - plot_lines_title .= "Posterior Density" - -} - --- Main -- - -main = do - rs <- runWithSystemRandom . generate $ conditionalLatentDistribution trns xy0 - - print ("Spike count: " ++ show (sum rs)) - - let rnbl = transducerRenderable rs - - void $ renderableToAspectWindow False 800 800 rnbl diff --git a/scripts/transducer.hs b/scripts/transducer.hs deleted file mode 100644 --- a/scripts/transducer.hs +++ /dev/null @@ -1,110 +0,0 @@ -{-# LANGUAGE TypeFamilies #-} - ---- Imports --- - - --- Goal -- - -import Goal.Core - -import Goal.Geometry -import Goal.Probability - - ---- Program --- - - --- Globals -- - -vr = 1 -mn = -4 -mx = 4 -nkrns = 10 -mus = range mn mx nkrns -sps = [ fromList Normal [mu,vr] | mu <- mus] - -gn1 = 2 -gn2 = 4 - -trns1 = modulateTransducerGain gn1 $ buildNormalTransducer sps -trns2 = modulateTransducerGain gn2 $ buildNormalTransducer sps - -x0 = 0 - -stps = 2000 -pltrng = range mn mx stps -laxprms = LinearAxisParams (show . round) 2 2 -iaxprms = LinearAxisParams show 3 3 -xaxprms = LinearAxisParams (show . round) 5 5 - - --- Functions -- - - --- Main -- - -main = do - - rs1 <- runWithSystemRandom . generate $ conditionalLatentDistribution trns1 x0 - rs2 <- runWithSystemRandom . generate $ conditionalLatentDistribution trns2 x0 - - let tclyt = execEC $ do - - layout_y_axis . laxis_generate .= scaledAxis laxprms (0,1.5) - layout_x_axis . laxis_generate .= autoScaledAxis xaxprms - --layout_y_axis . laxis_title .= "Activation" - layout_y_axis . laxis_override .= axisGridHide - - --layout_x_axis . laxis_title .= "Stimulus" - layout_x_axis . laxis_override .= axisGridHide - - plot . liftEC $ do - --plot_lines_title .= "Tuning Curves" - plot_lines_style .= solidLine 2 (opaque black) - plot_lines_values .= ( zip pltrng <$> transpose - (listCoordinates . (gn1 />) . potentialMapping <$> conditionalLatentDistributions trns1 pltrng) ) - - let rsplytfun trns rs = execEC $ do - - let posterior = conditionalObservableDistribution trns rs - scl = 10 - - --layoutlr_title .= ("μ=" ++ showFFloat (Just 3) mu "" ++ "; σ=" ++ showFFloat (Just 3) sd "") - - layoutlr_left_axis . laxis_generate .= scaledAxis laxprms (0,2) - --layoutlr_left_axis . laxis_title .= "Probability Density" - layoutlr_left_axis . laxis_override .= axisGridHide - - layoutlr_right_axis . laxis_generate .= scaledIntAxis iaxprms (0,round scl) - --layoutlr_right_axis . laxis_title .= "Response Count" - layoutlr_right_axis . laxis_override .= axisGridHide - - --layoutlr_x_axis . laxis_title .= "Stimulus" - layoutlr_margin .= 10 - - layoutlr_x_axis . laxis_override .= axisGridHide - layoutlr_x_axis . laxis_generate .= autoScaledAxis xaxprms - - layoutlr_plots - .= [ Left $ vlinePlot "" (solidLine 2 $ opaque black) x0 ] - - plotRight . liftEC $ do - plot_points_style .= filledCircles 3 (opaque black) - plot_points_values .= zip mus rs - --plot_points_title .= "Response" - - plotLeft . liftEC $ do - plot_lines_style .= solidLine 2 (opaque red) - plot_lines_values .= [zip pltrng $ density posterior <$> pltrng] - --plot_lines_title .= "Posterior Density" - - let rsplyt1 = rsplytfun trns1 rs1 - rsplyt2 = rsplytfun trns2 rs2 - rsplyt3 = rsplytfun trns2 (zipWith (+) rs1 rs2) - - let rnbl = toRenderable . weights (1,1) - $ tval (StackedLayouts [StackedLayout tclyt, StackedLayoutLR rsplyt2] True) - .|. tval (StackedLayouts [StackedLayoutLR rsplyt1, StackedLayoutLR rsplyt3] True) - - void $ renderableToAspectWindow False 1200 800 rnbl - --void $ renderableToFile (FileOptions (600,300) PDF) "population-code.pdf" rnbl diff --git a/scripts/univariate.hs b/scripts/univariate.hs deleted file mode 100644 --- a/scripts/univariate.hs +++ /dev/null @@ -1,99 +0,0 @@ -{-# LANGUAGE TypeOperators, TypeFamilies, FlexibleContexts #-} - ---- Imports --- - - --- Goal -- - -import Goal.Core -import Goal.Geometry -import Goal.Probability - - ---- Globals --- - -nsmps = 20 - --- Bernoulli -- - -(mnB,mxB) = (0,1) -bnsB = 2 -truB = chart Standard $ fromList Bernoulli [0.7] -toDoubleB = coordinate 0 . sufficientStatistic Bernoulli -rngB = [False,True] - --- Categorical -- - -(mnC,mxC) = (0,4) -bnsC = 5 -toDoubleC = fromIntegral -truC = chart Standard $ fromList (Categorical [0,1,2,3,4]) [0.1,0.4,0.1,0.2] -rngC = [0..4] - --- Poisson -- - -(mnP,mxP) = (0,20) -bnsP = 20 -toDoubleP = fromIntegral -truP = chart Standard $ fromList Poisson [5] -rngP = [0..20] - --- Normal -- - -(mnN,mxN) = (-3,7) -bnsN = 20 -toDoubleN = id -truN = chart Standard $ fromList Normal [2,0.7] -rngN = [-3,-2.99..7] - --- Layout -- - -generateLayout :: ( Show m, Transition Standard Mixture m, Transition Standard Natural m - , MaximumLikelihood Standard m, AbsolutelyContinuous Standard m, Generative Standard m , ExponentialFamily m ) - => Int -> Double -> Double -> (Sample m -> Double) -> [Sample m] -> Standard :#: m -> IO (LayoutLR Double Int Double) -generateLayout nb mn mx toDouble rng p = do - - let m = manifold p - lineFun1 p' = zip (toDouble <$> rng) $ density p' <$> rng - lineFun2 p' = zip (toDouble <$> rng) $ density p' <$> rng - - smps <- runWithSystemRandom . replicateM nsmps $ generate p - - let mle1 = chart Standard $ mle m smps - let hstplt = histogramPlot nb mn mx [toDouble <$> smps] . execEC $ do - plot_bars_titles .= ["Samples"] - plot_bars_item_styles .= [(solidFillStyle $ opaque blue, Nothing)] - - return . histogramLayoutLR hstplt . execEC $ do - - layoutlr_title .= (show (manifold p) ++ "; KLD: " ++ take 5 (showFFloat (Just 3) (klDivergence mle1 p) "")) - layoutlr_left_axis . laxis_title .= "Sample Count" - layoutlr_right_axis . laxis_title .= "Probability Mass" - layoutlr_x_axis . laxis_title .= "Value" - - plotRight . liftEC $ do - plot_lines_style .= dashedLine 3 [2,1] (opaque black) - plot_lines_title .= "True" - plot_lines_values .= [lineFun1 p] - - plotRight . liftEC $ do - plot_lines_style .= dashedLine 3 [10,5] (opaque red) - plot_lines_title .= "Standard MLE" - plot_lines_values .= [ lineFun1 mle1 ] - - plotRight . liftEC $ do - plot_lines_style .= dashedLine 3 [7,3] (opaque purple) - plot_lines_title .= "Exponential Family MLE" - plot_lines_values .= [ lineFun2 . chart Natural $ mle m smps ] - -main = do - - lytB <- tval <$> generateLayout bnsB mnB mxB toDoubleB rngB truB - lytC <- tval <$> generateLayout bnsC mnC mxC toDoubleC rngC truC - lytP <- tval <$> generateLayout bnsP mnP mxP toDoubleP rngP truP - lytN <- tval <$> generateLayout bnsN mnN mxN toDoubleN rngN truN - - let grd1 = lytB .|. lytC - grd2 = lytP .|. lytN - - renderableToAspectWindow False 800 600 . toRenderable . weights (1,1) $ grd1 ./. grd2