diff --git a/Goal/Probability.hs b/Goal/Probability.hs new file mode 100644 --- /dev/null +++ b/Goal/Probability.hs @@ -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 diff --git a/Goal/Probability/Distributions.hs b/Goal/Probability/Distributions.hs new file mode 100644 --- /dev/null +++ b/Goal/Probability/Distributions.hs @@ -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 + +-} diff --git a/Goal/Probability/ExponentialFamily.hs b/Goal/Probability/ExponentialFamily.hs new file mode 100644 --- /dev/null +++ b/Goal/Probability/ExponentialFamily.hs @@ -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 diff --git a/Goal/Probability/Graphical.hs b/Goal/Probability/Graphical.hs new file mode 100644 --- /dev/null +++ b/Goal/Probability/Graphical.hs @@ -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 + diff --git a/Goal/Probability/Graphical/Harmonium.hs b/Goal/Probability/Graphical/Harmonium.hs new file mode 100644 --- /dev/null +++ b/Goal/Probability/Graphical/Harmonium.hs @@ -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) diff --git a/Goal/Probability/Graphical/NeuralNetwork.hs b/Goal/Probability/Graphical/NeuralNetwork.hs new file mode 100644 --- /dev/null +++ b/Goal/Probability/Graphical/NeuralNetwork.hs @@ -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 + -} diff --git a/Goal/Probability/Statistical.hs b/Goal/Probability/Statistical.hs new file mode 100644 --- /dev/null +++ b/Goal/Probability/Statistical.hs @@ -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) + +-} diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/goal-probability.cabal b/goal-probability.cabal new file mode 100644 --- /dev/null +++ b/goal-probability.cabal @@ -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 diff --git a/scripts/backpropagation.hs b/scripts/backpropagation.hs new file mode 100644 --- /dev/null +++ b/scripts/backpropagation.hs @@ -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 + -} diff --git a/scripts/cross-entropy-descent.hs b/scripts/cross-entropy-descent.hs new file mode 100644 --- /dev/null +++ b/scripts/cross-entropy-descent.hs @@ -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 diff --git a/scripts/divergence.hs b/scripts/divergence.hs new file mode 100644 --- /dev/null +++ b/scripts/divergence.hs @@ -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 + + diff --git a/scripts/multivariate.hs b/scripts/multivariate.hs new file mode 100644 --- /dev/null +++ b/scripts/multivariate.hs @@ -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 diff --git a/scripts/poisson-binomial.hs b/scripts/poisson-binomial.hs new file mode 100644 --- /dev/null +++ b/scripts/poisson-binomial.hs @@ -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 diff --git a/scripts/transducer-field.hs b/scripts/transducer-field.hs new file mode 100644 --- /dev/null +++ b/scripts/transducer-field.hs @@ -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 diff --git a/scripts/transducer.hs b/scripts/transducer.hs new file mode 100644 --- /dev/null +++ b/scripts/transducer.hs @@ -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 diff --git a/scripts/univariate.hs b/scripts/univariate.hs new file mode 100644 --- /dev/null +++ b/scripts/univariate.hs @@ -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