goal-probability (empty) → 0.1
raw patch · 18 files changed
+2339/−0 lines, 18 filesdep +basedep +goal-coredep +goal-geometrysetup-changed
Dependencies added: base, goal-core, goal-geometry, goal-probability, hmatrix, math-functions, mwc-random, mwc-random-monad, statistics, vector
Files
- Goal/Probability.hs +73/−0
- Goal/Probability/Distributions.hs +650/−0
- Goal/Probability/ExponentialFamily.hs +98/−0
- Goal/Probability/Graphical.hs +9/−0
- Goal/Probability/Graphical/Harmonium.hs +214/−0
- Goal/Probability/Graphical/NeuralNetwork.hs +239/−0
- Goal/Probability/Statistical.hs +131/−0
- LICENSE +30/−0
- Setup.hs +2/−0
- goal-probability.cabal +134/−0
- scripts/backpropagation.hs +120/−0
- scripts/cross-entropy-descent.hs +114/−0
- scripts/divergence.hs +81/−0
- scripts/multivariate.hs +100/−0
- scripts/poisson-binomial.hs +51/−0
- scripts/transducer-field.hs +84/−0
- scripts/transducer.hs +110/−0
- scripts/univariate.hs +99/−0
+ Goal/Probability.hs view
@@ -0,0 +1,73 @@+module Goal.Probability+ ( module System.Random.MWC+ , module System.Random.MWC.Monad+ , module Goal.Probability.Statistical+ , module Goal.Probability.ExponentialFamily+ , module Goal.Probability.Distributions+ , module Goal.Probability.Graphical+ , module Goal.Probability.Graphical.Harmonium+ , module Goal.Probability.Graphical.NeuralNetwork+ , module Goal.Probability+ ) where+++--- Imports ---+++-- 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.Distributions+import Goal.Probability.Graphical+import Goal.Probability.Graphical.Harmonium+import Goal.Probability.Graphical.NeuralNetwork++-- Package --++import Goal.Core+import Goal.Geometry++--- Stochastic Functions ---+++seed :: RandST s Seed+-- | This little guy creates a seed. It's necessary to avoid name space+-- collisions.+seed = S.save++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++noisyFunction :: (Generative c m, Num (Sample m))+ => (c :#: m) -- ^ Noise model+ -> (x -> Sample m) -- ^ Function+ -> x+ -> RandST r (Sample m)+-- | Returns a sample from the given function with added noise.+noisyFunction m f x = do+ ns <- generate 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
+ Goal/Probability/Distributions.hs view
@@ -0,0 +1,650 @@+-- | 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+ ) where++-- Package --++import Goal.Core+import Goal.Probability.Statistical+import Goal.Probability.ExponentialFamily++import Goal.Geometry++-- Qualified --++import qualified Data.Vector.Storable as C+import qualified Numeric.LinearAlgebra.HMatrix as M++-- Unqualified --++import System.Random.MWC.Monad+import System.Random.MWC.Distributions.Monad+import Statistics.Sample hiding (mean)+import Numeric.SpecFunctions++-- Uniform --++data Uniform = Uniform Double Double deriving (Eq, Read, Show)++instance Manifold Uniform where+ dimension _ = 0++instance Statistical Uniform where+ type SampleSpace Uniform = Continuum+ sampleSpace _ = Continuum++instance Generative Standard Uniform where+ generate p =+ let (Uniform a b) = manifold p+ in uniformR (a,b)++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++-- Bernoulli Distribution --++-- | The Bernoulli 'Family' with 'SampleSpace' 'Bernoulli' = 'Bool' (because why not).+data Bernoulli = Bernoulli deriving (Eq, Read, Show)++instance Manifold Bernoulli where+ dimension _ = 1++instance Statistical Bernoulli where+ type SampleSpace Bernoulli = Boolean+ sampleSpace Bernoulli = Boolean++instance Generative Standard Bernoulli where+ generate p = bernoulli . C.head $ coordinates p++instance AbsolutelyContinuous Standard Bernoulli where+ density p True = C.head $ coordinates p+ density p False = 1 - C.head (coordinates p)++instance MaximumLikelihood Standard Bernoulli where+ mle _ bls = fromList Bernoulli [mean $ toDouble <$> bls]+ where toDouble True = 1+ toDouble False = 0++instance Legendre Natural Bernoulli where+ potential p = log $ 1 + exp (coordinate 0 p)+ potentialDifferentials p = fromList (Tangent p) [logistic $ coordinate 0 p]++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]++instance ExponentialFamily Bernoulli where+ baseMeasure _ _ = 1+ sufficientStatistic Bernoulli True = fromList Bernoulli [1]+ sufficientStatistic Bernoulli False = fromList Bernoulli [0]++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)]++instance Transition Standard Mixture Bernoulli where+ transition = breakChart++instance Transition Mixture Standard Bernoulli where+ transition = breakChart++instance Transition Standard Natural Bernoulli where+ transition = potentialMapping . chart Mixture . transition++instance Transition Natural Standard Bernoulli where+ transition = transition . potentialMapping++instance Generative Natural Bernoulli where+ generate = standardGenerate+++-- Binomial Distribution --++newtype Binomial = Binomial { binomialTrials :: Int } deriving (Eq, Read, Show)++instance Manifold Binomial where+ dimension _ = 1++instance Statistical Binomial where+ type SampleSpace Binomial = [Int]+ sampleSpace (Binomial n) = [0..n]++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 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 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)]+++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 ExponentialFamily Binomial where+ baseMeasure (Binomial n) = choose n+ sufficientStatistic s k = fromList s [fromIntegral k]++instance Transition Standard Natural Binomial where+ transition = potentialMapping . chart Mixture . transition++instance Transition Natural Standard Binomial where+ transition = chart Standard . transition . potentialMapping++instance Transition Standard Mixture Binomial where+ transition p = breakChart $ alterCoordinates (* (fromIntegral . binomialTrials $ manifold p)) p++instance Transition Mixture Standard Binomial where+ transition p = breakChart $ alterCoordinates (/ (fromIntegral . binomialTrials $ manifold p)) p++-- Categorical Distribution --++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.++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 Discrete s => Manifold (Categorical s) where+ dimension (Categorical s) = length (elements s) - 1++instance Discrete s => Statistical (Categorical s) where+ type SampleSpace (Categorical s) = s+ sampleSpace (Categorical ks) = ks++instance Discrete s => Generative Standard (Categorical s) where+ generate p = generateCategorical (samples $ manifold p) (coordinates 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 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 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 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 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 Discrete s => Transition Standard Mixture (Categorical s) where+ transition = breakChart++instance Discrete s => Transition Mixture Standard (Categorical s) where+ transition = breakChart++instance Discrete s => Transition Standard Natural (Categorical s) where+ transition = potentialMapping . chart Mixture . transition++instance Discrete s => Transition Natural Standard (Categorical s) where+ transition = transition . potentialMapping++-- Curved Categorical Distribution --++newtype CurvedCategorical s = CurvedCategorical s deriving (Show,Eq,Read)++instance Discrete s => Manifold (CurvedCategorical s) where+ dimension = length . samples++instance Discrete s => Statistical (CurvedCategorical s) where+ type SampleSpace (CurvedCategorical s) = s+ sampleSpace (CurvedCategorical s) = s++instance Discrete s => Generative Standard (CurvedCategorical s) where+ generate p = generateCategorical (samples $ manifold p) (coordinates p)++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++-- Poisson Distribution --++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)++data Poisson = Poisson deriving (Eq, Read, Show)++instance Manifold Poisson where+ dimension _ = 1++instance Statistical Poisson where+ type SampleSpace Poisson = NaturalNumbers+ sampleSpace _ = NaturalNumbers++instance Generative Standard Poisson where+ generate d = generatePoisson . C.head $ coordinates d++instance AbsolutelyContinuous Standard Poisson where+ density d k =+ let ps = coordinates d+ lmda = C.head ps+ in lmda^k / factorial k * exp (-lmda)++instance MaximumLikelihood Standard Poisson where+ mle _ xs = fromList Poisson . (:[]) . mean $ fromIntegral <$> xs++instance ExponentialFamily Poisson where+ sufficientStatistic Poisson = fromCoordinates Poisson . C.singleton . fromIntegral+ baseMeasure _ k = recip $ factorial k++instance Legendre Natural Poisson where+ potential p = exp $ coordinate 0 p+ potentialDifferentials p = fromCoordinates (Tangent p) . exp $ coordinates p++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]++instance Transition Standard Natural Poisson where+ transition = transition . chart Mixture . transition++instance Transition Natural Standard Poisson where+ transition = transition . potentialMapping++instance Transition Standard Mixture Poisson where+ transition = breakChart++instance Transition Mixture Standard Poisson where+ transition = breakChart++instance Generative Natural Poisson where+ generate = standardGenerate++-- Normal Distribution --++data Normal = Normal deriving (Show,Eq,Read)++instance Manifold Normal where+ dimension _ = 2++instance Statistical Normal where+ type SampleSpace Normal = Continuum+ sampleSpace _ = Continuum++instance Generative Standard Normal where+ generate p =+ let [mu,vr] = listCoordinates p+ in normal mu $ sqrt vr++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 MaximumLikelihood Standard Normal where+ mle _ xs =+ let (mu,vr) = meanVariance $ C.fromList xs+ in fromList Normal [mu,vr]++instance ExponentialFamily Normal where+ sufficientStatistic Normal x = fromList Normal [x,x**2]+ baseMeasure _ _ = recip . sqrt $ 2 * pi++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 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 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 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 Transition Mixture Standard Normal where+ transition p =+ let [eta0,eta1] = listCoordinates p+ in fromList Normal [eta0, eta1 - eta0^2]++instance Transition Standard Natural Normal where+ transition p =+ let [mu,vr] = listCoordinates p+ in fromList Normal [mu / vr, negate . recip $ 2 * vr]++instance Transition Natural Standard Normal where+ transition p =+ let [tht0,tht1] = listCoordinates p+ in fromList Normal [-0.5 * tht0 / tht1, negate . recip $ 2 * tht1]++instance Generative Natural Normal where+ generate = standardGenerate++-- MeanNormal Distribution --++data MeanNormal = MeanNormal Double deriving (Show,Eq,Read)++instance Manifold MeanNormal where+ dimension _ = 1+++instance Statistical MeanNormal where+ type SampleSpace MeanNormal = Continuum+ sampleSpace _ = Continuum++instance Generative Standard MeanNormal where+ generate p = do+ let (MeanNormal vr) = manifold p+ normal (coordinate 0 p) $ sqrt vr++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 MaximumLikelihood Standard MeanNormal where+ mle mnrm xs = fromList mnrm [mean xs]++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 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 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 Riemannian Natural MeanNormal where+ metric p =+ let (MeanNormal vr) = manifold p+ in fromList (Tensor (Tangent p) (Tangent p)) [vr]++instance Transition Standard Natural MeanNormal where+ transition = potentialMapping . chart Mixture . breakChart++instance Transition Natural Standard MeanNormal where+ transition = breakChart . potentialMapping++instance Transition Standard Mixture MeanNormal where+ transition = breakChart++instance Transition Mixture Standard MeanNormal where+ transition = breakChart++-- Multivariate Normal --++data MultivariateNormal = MultivariateNormal { sampleSpaceDimension :: Int } deriving (Eq, Read, Show)++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++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++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 Manifold MultivariateNormal where+ dimension (MultivariateNormal n) = n + n^2++instance Statistical MultivariateNormal where+ type SampleSpace MultivariateNormal = Euclidean+ sampleSpace (MultivariateNormal n) = Euclidean n++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 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 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 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 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 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 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 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 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 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)+++{-+--- Graveyard ---+++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++-- Exponential Distribution --++data Exponential = Exponential deriving (Eq,Read,Show)++instance Manifold Exponential where+ dimension _ = 1++type instance SampleSpace Exponential = Continuum++instance Statistical Exponential where+ sampleSpace _ = Continuum++instance Generative Standard Exponential where+ generate = exponential . C.head . coordinates++instance AbsolutelyContinuous Standard Exponential where+ density p x =+ let lmda = C.head $ coordinates p+ in lmda * exp (negate $ lmda * x)++instance MaximumLikelihood Standard Exponential where+ mle _ xs = chart Standard . fromList Exponential . (:[]) . recip . mean $ xs++instance Legendre Natural Exponential where+ potential p = negate . log . negate $ coordinate 0 p+ potentialDifferentials p = fromCoordinates (Tangent p) . negate $ coordinates p++instance Legendre Mixture Exponential where+ potential p = 1 - log eta+ potentialDifferentials p =++instance ExponentialFamily Exponential where+ sufficientStatistic Exponential = fromCoordinates Exponential . C.singleton+ baseMeasure _ _ = 1++instance Transition Standard Natural Exponential where+ transition = breakChart . alterCoordinates negate++instance Transition Natural Standard Exponential where+ transition = breakChart . alterCoordinates negate++-}
+ Goal/Probability/ExponentialFamily.hs view
@@ -0,0 +1,98 @@+module Goal.Probability.ExponentialFamily (+ -- * Exponential Families+ ExponentialFamily (sufficientStatistic, baseMeasure)+ , sufficientStatisticN+ -- ** Dual Parameters+ , Natural (Natural)+ , Mixture (Mixture)+ -- ** Divergence+ , klDivergence+ , relativeEntropy+ ) where++--- Imports ---+++-- Package --++import Goal.Probability.Statistical++import Goal.Geometry+++--- 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.+--+-- '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++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)++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)++relativeEntropy+ :: (ExponentialFamily m, Transition c Mixture m, Transition d Natural m)+ => c :#: m -> d :#: m -> Double+relativeEntropy p q = klDivergence q p++-- | A parameterization in terms of the natural coordinates of an exponential family.+data Natural = Natural++-- | A representation in terms of the mean sufficient statistics of an exponential family.+data Mixture = Mixture++instance Primal Natural where+ type Dual Natural = Mixture++instance Primal Mixture where+ type Dual Mixture = Natural+++--- Instances ---+++-- Generic --++instance ExponentialFamily m => MaximumLikelihood Mixture m where+ mle = sufficientStatisticN++instance ExponentialFamily m => MaximumLikelihood Natural m where+ mle m xs = potentialMapping $ sufficientStatisticN m xs++-- Replicated --++instance ExponentialFamily m => ExponentialFamily (Replicated m) where+ sufficientStatistic (Replicated m _) xs =+ joinReplicated $ sufficientStatistic m <$> xs+ baseMeasure (Replicated m _) xs = product $ baseMeasure m <$> xs++-- Fisher Manifolds --++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 ExponentialFamily m => Transition Mixture Natural m where+ transition = potentialMapping++instance ExponentialFamily m => Transition Natural Mixture m where+ transition = potentialMapping
+ Goal/Probability/Graphical.hs view
@@ -0,0 +1,9 @@+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+
+ Goal/Probability/Graphical/Harmonium.hs view
@@ -0,0 +1,214 @@+-- | 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)
+ Goal/Probability/Graphical/NeuralNetwork.hs view
@@ -0,0 +1,239 @@+-- | 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+ -}
+ Goal/Probability/Statistical.hs view
@@ -0,0 +1,131 @@+module Goal.Probability.Statistical (+ -- * Stastical Manifolds+ Statistical (sampleSpace)+ , Sample+ , samples+ , SampleSpace+ -- ** Standard Chart+ , Standard (Standard)+ , standardGenerate+ -- ** Distributions+ , Generative (generate)+ , AbsolutelyContinuous (density)+ , expectation+ , MaximumLikelihood (mle)+ ) where+++--- Imports ---+++-- Package --++import Goal.Geometry++-- Unqualified --++import System.Random.MWC.Monad+++--- Test Bed ---+++--- 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 'Sample' is an 'Element' of the 'SampleSpace'.+type Sample m = Element (SampleSpace m)++samples :: (Discrete (SampleSpace m), Statistical m) => m -> [Sample m]+-- | The list of 'Sample's.+samples = elements . sampleSpace++-- | 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)++-- | 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++-- | '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+expectation p f =+ let xs = elements . sampleSpace $ manifold p+ in sum $ zipWith (*) (f <$> xs) (density p <$> xs)+++-- | 'mle' computes the 'MaximumLikelihood' estimator.+class Statistical m => MaximumLikelihood c m where+ mle :: m -> [Sample m] -> c :#: m++-- Standard Chart --++-- | 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++--- Instances ---+++-- DirectSums --++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 (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 (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++-- 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 m, Generative c m) => Generative c (Replicated m) where+ generate = sequence . mapReplicated generate++instance (Statistical m, AbsolutelyContinuous Standard m) => AbsolutelyContinuous Standard (Replicated m) where+ density ds xs = product $ zipWith ($) (mapReplicated density ds) xs++instance (Statistical m, Transition Standard c m) => Transition Standard c (Replicated m) where+ transition = joinReplicated . mapReplicated transition++instance (Statistical m, Transition c Standard m) => Transition c Standard (Replicated m) where+ transition = joinReplicated . mapReplicated transition+++--- Graveyard ---+++{-+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)++-}
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c) 2014, Sacha Sokoloski++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++ * Redistributions in binary form must reproduce the above+ copyright notice, this list of conditions and the following+ disclaimer in the documentation and/or other materials provided+ with the distribution.++ * Neither the name of Sacha Sokoloski nor the names of other+ contributors may be used to endorse or promote products derived+ from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ goal-probability.cabal view
@@ -0,0 +1,134 @@+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+license-file: LICENSE+author: Sacha Sokoloski+maintainer: sokolo@mis.mpg.de+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+ build-depends:+ base==4.*,+ goal-core==0.1,+ goal-geometry==0.1,+ goal-probability==0.1+ default-language: Haskell2010++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+ build-depends:+ base==4.*,+ goal-core==0.1,+ goal-geometry==0.1,+ goal-probability==0.1+ default-language: Haskell2010++executable backpropagation+ 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+ build-depends:+ base==4.*,+ goal-core==0.1,+ goal-geometry==0.1,+ goal-probability==0.1+ default-language: Haskell2010
+ scripts/backpropagation.hs view
@@ -0,0 +1,120 @@+{-# 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+ -}
+ scripts/cross-entropy-descent.hs view
@@ -0,0 +1,114 @@+{-# 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
+ scripts/divergence.hs view
@@ -0,0 +1,81 @@+{-# 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++
+ scripts/multivariate.hs view
@@ -0,0 +1,100 @@+--- 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
+ scripts/poisson-binomial.hs view
@@ -0,0 +1,51 @@+-- 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
+ scripts/transducer-field.hs view
@@ -0,0 +1,84 @@+--- 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
+ scripts/transducer.hs view
@@ -0,0 +1,110 @@+{-# 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
+ scripts/univariate.hs view
@@ -0,0 +1,99 @@+{-# 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