goal-probability 0.1 → 0.20

raw patch · 22 files changed

+2337/−2015 lines, 22 filesdep +bytestringdep +cassavadep +containersdep −math-functionsdep −mwc-random-monaddep ~basedep ~goal-coredep ~goal-geometry

Dependencies added: bytestring, cassava, containers, criterion, ghc-typelits-knownnat, ghc-typelits-natnormalise, hmatrix-special, parallel

Dependencies removed: math-functions, mwc-random-monad

Dependency ranges changed: base, goal-core, goal-geometry, goal-probability, hmatrix, mwc-random, statistics, vector

Files

• Goal/Probability.hs +204/−44
• Goal/Probability/Conditional.hs +228/−0
• Goal/Probability/Distributions.hs +456/−493
• Goal/Probability/Distributions/CoMPoisson.hs +195/−0
• Goal/Probability/Distributions/Gaussian.hs +560/−0
• Goal/Probability/ExponentialFamily.hs +194/−61
• Goal/Probability/Graphical.hs +0/−9
• Goal/Probability/Graphical/Harmonium.hs +0/−214
• Goal/Probability/Graphical/NeuralNetwork.hs +0/−239
• Goal/Probability/Statistical.hs +173/−79
• README.md +72/−0
• benchmarks/backpropagation.hs +79/−0
• benchmarks/regression.hs +112/−0
• goal-probability.cabal +64/−117
• scripts/backpropagation.hs +0/−120
• scripts/cross-entropy-descent.hs +0/−114
• scripts/divergence.hs +0/−81
• scripts/multivariate.hs +0/−100
• scripts/poisson-binomial.hs +0/−51
• scripts/transducer-field.hs +0/−84
• scripts/transducer.hs +0/−110
• scripts/univariate.hs +0/−99

@@ -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)

@@ -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]+++

@@ -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

@@ -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).+-- (<https://rss.onlinelibrary.wiley.com/doi/abs/10.1111/j.1467-9876.2005.00474.x>) 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

@@ -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++

@@ -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

@@ -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-

@@ -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)

@@ -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-    -}

@@ -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

@@ -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).+

@@ -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 ]

@@ -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

@@ -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

@@ -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-    -}

@@ -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

@@ -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--

@@ -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

@@ -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

@@ -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

@@ -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

@@ -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