diff --git a/Goal/Graphical.hs b/Goal/Graphical.hs
new file mode 100644
--- /dev/null
+++ b/Goal/Graphical.hs
@@ -0,0 +1,31 @@
+{-# OPTIONS_GHC -fplugin=GHC.TypeLits.KnownNat.Solver -fplugin=GHC.TypeLits.Normalise -fconstraint-solver-iterations=10 #-}
+{-# LANGUAGE
+    RankNTypes,
+    TypeOperators,
+    FlexibleContexts,
+    ScopedTypeVariables
+#-}
+-- | The main module of goal-probability. Import this module to use all the
+-- types, functions, and classes provided by goal-probability.
+module Goal.Graphical
+    ( -- * Package Exports
+      module Goal.Graphical.Models
+    , module Goal.Graphical.Models.Dynamic
+    , module Goal.Graphical.Models.Harmonium
+    , module Goal.Graphical.Models.Harmonium.FactorAnalysis
+    , module Goal.Graphical.Learning
+    , module Goal.Graphical.Inference
+    ) where
+
+
+--- Imports ---
+
+
+-- Re-exports --
+
+import Goal.Graphical.Models
+import Goal.Graphical.Models.Dynamic
+import Goal.Graphical.Models.Harmonium
+import Goal.Graphical.Models.Harmonium.FactorAnalysis
+import Goal.Graphical.Learning
+import Goal.Graphical.Inference
diff --git a/Goal/Graphical/Inference.hs b/Goal/Graphical/Inference.hs
new file mode 100644
--- /dev/null
+++ b/Goal/Graphical/Inference.hs
@@ -0,0 +1,169 @@
+{-# OPTIONS_GHC -fplugin=GHC.TypeLits.KnownNat.Solver -fplugin=GHC.TypeLits.Normalise -fconstraint-solver-iterations=10 #-}
+{-# LANGUAGE
+    RankNTypes,
+    PolyKinds,
+    DataKinds,
+    TypeOperators,
+    FlexibleContexts,
+    FlexibleInstances,
+    TypeApplications,
+    ScopedTypeVariables,
+    TypeFamilies
+#-}
+-- | Infering latent variables in graphical models.
+module Goal.Graphical.Inference
+    ( -- * Inference
+      conjugatedBayesRule
+    -- * Recursive
+    , conjugatedRecursiveBayesianInference
+    -- * Dynamic
+    , conjugatedPredictionStep
+    , conjugatedForwardStep
+    -- * Conjugation
+    , regressConjugationParameters
+    , conjugationCurve
+    ) where
+
+--- Imports ---
+
+
+-- Goal --
+
+import Goal.Core
+import Goal.Geometry
+import Goal.Probability
+
+import Goal.Graphical.Models.Harmonium
+
+import qualified Goal.Core.Vector.Storable as S
+
+import Data.List
+
+
+--- Inference ---
+
+
+-- | The posterior distribution given a prior and likelihood, where the
+-- likelihood is conjugated.
+conjugatedBayesRule
+    :: forall f y x z w
+    . ( Map Natural f x y, Bilinear f y x, ConjugatedLikelihood f y x z w )
+    => Natural # Affine f y z x
+    -> Natural # w
+    -> SamplePoint z
+    -> Natural # w
+conjugatedBayesRule lkl prr z =
+    let pstr = fst . split . transposeHarmonium $ joinConjugatedHarmonium lkl prr
+        mz :: Mean # z
+        mz = sufficientStatistic z
+     in pstr >.+> mz
+
+
+--- Recursive ---
+
+
+-- | The posterior distribution given a prior and likelihood, where the
+-- likelihood is conjugated.
+conjugatedRecursiveBayesianInference
+    :: ( Map Natural f x y, Bilinear f y x, ConjugatedLikelihood f y x z w )
+    => Natural # Affine f y z x -- ^ Likelihood
+    -> Natural # w -- ^ Prior
+    -> Sample z -- ^ Observations
+    -> [Natural # w] -- ^ Updated prior
+conjugatedRecursiveBayesianInference lkl = scanl' (conjugatedBayesRule lkl)
+
+
+-- Dynamical ---
+
+
+-- | The predicted distribution given a current distribution and transition
+-- distribution, where the transition distribution is (doubly) conjugated.
+conjugatedPredictionStep
+    :: (ConjugatedLikelihood f x x w w, Bilinear f x x)
+    => Natural # Affine f x w x
+    -> Natural # w
+    -> Natural # w
+conjugatedPredictionStep trns prr =
+    snd . splitConjugatedHarmonium . transposeHarmonium
+        $ joinConjugatedHarmonium trns prr
+
+-- | Forward inference based on conjugated models: priors at a previous time are
+-- first predicted into the current time, and then updated with Bayes rule.
+conjugatedForwardStep
+    :: ( ConjugatedLikelihood g x x w w, Bilinear g x x
+       , ConjugatedLikelihood f y x z w, Bilinear f y x
+       , Map Natural g x x, Map Natural f x y )
+    => Natural # Affine g x w x -- ^ Transition Distribution
+    -> Natural # Affine f y z x -- ^ Emission Distribution
+    -> Natural # w -- ^ Beliefs at time $t-1$
+    -> SamplePoint z -- ^ Observation at time $t$
+    -> Natural # w -- ^ Beliefs at time $t$
+conjugatedForwardStep trns emsn prr z =
+    flip (conjugatedBayesRule emsn) z $ conjugatedPredictionStep trns prr
+
+
+--- Approximate Conjugation ---
+
+
+-- | Computes the conjugation curve given a set of conjugation parameters,
+-- at the given set of points.
+conjugationCurve
+    :: ExponentialFamily x
+    => Double -- ^ Conjugation shift
+    -> Natural # x -- ^ Conjugation parameters
+    -> Sample x -- ^ Samples points
+    -> [Double] -- ^ Conjugation curve at sample points
+conjugationCurve rho0 rprms mus = (\x -> rprms <.> sufficientStatistic x + rho0) <$> mus
+
+-- Linear Least Squares
+
+-- | Returns the conjugation parameters which best satisfy the conjugation
+-- equation for the given population code according to linear regression.
+regressConjugationParameters
+    :: (Map Natural f z x, LegendreExponentialFamily z, ExponentialFamily x)
+    => Natural # f z x -- ^ PPC
+    -> Sample x -- ^ Sample points
+    -> (Double, Natural # x) -- ^ Approximate conjugation parameters
+regressConjugationParameters lkl mus =
+    let dpnds = potential <$> lkl >$>* mus
+        indpnds = independentVariables0 lkl mus
+        (rho0,rprms) = S.splitAt $ S.linearLeastSquares indpnds dpnds
+     in (S.head rho0, Point rprms)
+
+--- Internal ---
+
+independentVariables0
+    :: forall f x z . ExponentialFamily x
+    => Natural # f z x
+    -> Sample x
+    -> [S.Vector (Dimension x + 1) Double]
+independentVariables0 _ mus =
+    let sss :: [Mean # x]
+        sss = sufficientStatistic <$> mus
+     in (S.singleton 1 S.++) . coordinates <$> sss
+
+
+---- | The posterior distribution given a prior and likelihood, where the
+---- posterior is normalized via numerical integration.
+--numericalRecursiveBayesianInference
+--    :: forall f z x .
+--        ( Map Natural f x z, Map Natural f z x, Bilinear f z x
+--        , LegendreExponentialFamily z, ExponentialFamily x, SamplePoint x ~ Double)
+--    => Double -- ^ Integral error bound
+--    -> Double -- ^ Sample space lower bound
+--    -> Double -- ^ Sample space upper bound
+--    -> Sample x -- ^ Centralization samples
+--    -> [Natural # Affine f z x] -- ^ Likelihoods
+--    -> Sample z -- ^ Observations
+--    -> (Double -> Double) -- ^ Prior
+--    -> (Double -> Double, Double) -- ^ Posterior Density and Log-Partition Function
+--numericalRecursiveBayesianInference errbnd mnx mxx xsmps lkls zs prr =
+--    let logbm = logBaseMeasure (Proxy @ x)
+--        logupst0 x lkl z =
+--            (z *<.< snd (splitAffine lkl)) <.> sufficientStatistic x - potential (lkl >.>* x)
+--        logupst x = sum $ logbm x : log (prr x) : zipWith (logupst0 x) lkls zs
+--        logprt = logIntegralExp errbnd logupst mnx mxx xsmps
+--        dns x = exp $ logupst x - logprt
+--     in (dns,logprt)
+
+
diff --git a/Goal/Graphical/Learning.hs b/Goal/Graphical/Learning.hs
new file mode 100644
--- /dev/null
+++ b/Goal/Graphical/Learning.hs
@@ -0,0 +1,201 @@
+{-# OPTIONS_GHC -fplugin=GHC.TypeLits.KnownNat.Solver -fplugin=GHC.TypeLits.Normalise -fconstraint-solver-iterations=10 #-}
+{-# LANGUAGE Arrows #-}
+-- | A collection of algorithms for optimizing harmoniums.
+
+module Goal.Graphical.Learning
+    ( -- * Expectation Maximization
+      expectationMaximization
+    , expectationMaximizationAscent
+    , gibbsExpectationMaximization
+    , latentProcessExpectationMaximization
+    , latentProcessExpectationMaximizationAscent
+    -- * Differentials
+    , harmoniumInformationProjectionDifferential
+    , contrastiveDivergence
+    ) where
+
+
+--- Imports ---
+
+
+-- Goal --
+
+import Goal.Core
+import Goal.Geometry
+import Goal.Probability
+
+import Goal.Graphical.Models
+import Goal.Graphical.Models.Harmonium
+import Goal.Graphical.Models.Dynamic
+
+
+--- Differentials ---
+
+
+-- | The differential of the dual relative entropy. Minimizing this results in
+-- the information projection of the model against the marginal distribution of
+-- the given harmonium. This is more efficient than the generic version.
+harmoniumInformationProjectionDifferential
+    :: ( Map Natural f y x, LegendreExponentialFamily z
+       , SamplePoint w ~ SamplePoint x, Translation z y
+       , ExponentialFamily x, ExponentialFamily w, Generative Natural w )
+    => Int
+    -> Natural # AffineHarmonium f y x z w -- ^ Harmonium
+    -> Natural # w -- ^ Model Distribution
+    -> Random (Mean # w) -- ^ Differential Estimate
+harmoniumInformationProjectionDifferential n hrm px = do
+    xs <- sample n px
+    let (lkl,nw) = split hrm
+        mys0 = lkl >$>* xs
+        mws = sufficientStatistic <$> xs
+        mys = zipWith (\mw my0 -> mw <.> (px - nw) - potential my0) mws mys0
+        ln = fromIntegral $ length xs
+        mwht = average mws
+        myht = sum mys / ln
+        foldfun (mw,my) (k,z0) = (k+1,z0 + ((my - myht) .> (mw - mwht)))
+    return . uncurry (/>) . foldr foldfun (-1,0) $ zip mws mys
+
+-- | Contrastive divergence on harmoniums (<https://www.mitpressjournals.org/doi/abs/10.1162/089976602760128018?casa_token=x_Twj1HaXcMAAAAA:7-Oq181aubCFwpG-f8Lo1wRKvGnmujzl8zjn9XbeO5nGhfvKCCQjsu4K4pJCkMNYUYWqc2qG7TRXBg Hinton, 2019>).
+contrastiveDivergence
+    :: ( Generative Natural z, ExponentialFamily z, Translation w x
+       , Generative Natural w, ExponentialFamily y, Translation z y
+       , LegendreExponentialFamily w, Bilinear f y x, Map Natural f x y
+       , Map Natural f y x, SamplePoint y ~ SamplePoint z
+       , SamplePoint x ~ SamplePoint w, ExponentialFamily x )
+      => Int -- ^ The number of contrastive divergence steps
+      -> Sample z -- ^ The initial states of the Gibbs chains
+      -> Natural # AffineHarmonium f y x z w -- ^ The harmonium
+      -> Random (Mean # AffineHarmonium f y x z w) -- ^ The gradient estimate
+contrastiveDivergence cdn zs hrm = do
+    xzs0 <- initialPass hrm zs
+    xzs1 <- iterateM' cdn (gibbsPass hrm) xzs0
+    return $ stochasticRelativeEntropyDifferential xzs0 xzs1
+
+
+--- Expectation Maximization ---
+
+
+-- | A single iteration of EM for 'Harmonium' based models.
+expectationMaximization
+    :: ( DuallyFlatExponentialFamily (AffineHarmonium f y x z w)
+       , ExponentialFamily z, Map Natural f x y, Bilinear f y x
+       , Translation z y, Translation w x, LegendreExponentialFamily w )
+    => Sample z
+    -> Natural # AffineHarmonium f y x z w
+    -> Natural # AffineHarmonium f y x z w
+expectationMaximization zs hrm = transition $ expectationStep zs hrm
+
+-- | Ascent of the EM objective on harmoniums for when the expectation
+-- step can't be computed in closed-form. The convergent harmonium distribution
+-- of the output harmonium-list is the result of 1 iteration of the EM
+-- algorithm.
+expectationMaximizationAscent
+    :: ( LegendreExponentialFamily (AffineHarmonium f y x z w)
+       , ExponentialFamily z, Map Natural f x y, Bilinear f y x
+       , Translation z y, Translation w x, LegendreExponentialFamily w )
+    => Double
+    -> GradientPursuit
+    -> Sample z
+    -> Natural # AffineHarmonium f y x z w
+    -> [Natural # AffineHarmonium f y x z w]
+expectationMaximizationAscent eps gp zs nhrm =
+    let mhrm' = expectationStep zs nhrm
+     in vanillaGradientSequence (relativeEntropyDifferential mhrm') (-eps) gp nhrm
+
+-- | Ascent of the EM objective on harmoniums for when the expectation
+-- step can't be computed in closed-form. The convergent harmonium distribution
+-- of the output harmonium-list is the result of 1 iteration of the EM
+-- algorithm.
+gibbsExpectationMaximization
+    :: ( ExponentialFamily z, Map Natural f x y, Manifold w, Map Natural f y x
+       , Translation z y, Translation w x, SamplePoint y ~ SamplePoint z
+       , SamplePoint w ~ SamplePoint x
+       , ExponentialFamily y, Generative Natural w, ExponentialFamily x
+       , Generative Natural z, Manifold (AffineHarmonium f y x z w)
+       , Bilinear f y x, LegendreExponentialFamily w )
+    => Double
+    -> Int
+    -> Int
+    -> GradientPursuit
+    -> Sample z -- ^ Observations
+    -> Natural # AffineHarmonium f y x z w -- ^ Current Harmonium
+    -> Chain Random (Natural # AffineHarmonium f y x z w) -- ^ Harmonium Chain
+gibbsExpectationMaximization eps cdn nbtch gp zs0 nhrm0 =
+    let mhrm0 = expectationStep zs0 nhrm0
+     in chainCircuit nhrm0 $ proc nhrm -> do
+         zs <- minibatcher nbtch zs0 -< ()
+         xzs0 <- arrM (uncurry initialPass) -< (nhrm,zs)
+         xzs1 <- arrM (\(x,y) -> iterateM' cdn (gibbsPass x) y) -< (nhrm,xzs0)
+         let dff = mhrm0 - averageSufficientStatistic xzs1
+         gradientCircuit eps gp -< (nhrm,vanillaGradient dff)
+
+latentProcessExpectationStep
+    :: ( ConjugatedLikelihood g x x w w, ConjugatedLikelihood f y x z w
+       , Transition Natural Mean w, Transition Natural Mean (AffineHarmonium g x x w w)
+       , Manifold (AffineHarmonium g x x w w)
+       , Bilinear g x x, Map Natural f x y, Bilinear f y x
+       , SamplePoint y ~ SamplePoint z )
+    => Observations (LatentProcess f g y x z w)
+    -> Natural # LatentProcess f g y x z w
+    -> (Mean # w, Mean # AffineHarmonium f y x z w, Mean # AffineHarmonium g x x w w)
+latentProcessExpectationStep zss ltnt =
+    let (prr,emsn,trns) = splitLatentProcess ltnt
+        (smthss,hrmss) = unzip $ conjugatedSmoothing0 prr emsn trns <$> zss
+        mprr = average $ toMean . head <$> smthss
+        mtrns = average $ toMean <$> concat hrmss
+        mws = toMean <$> concat smthss
+        mzs = sufficientStatistic <$> concat zss
+        mys = anchor <$> mzs
+        mxs = anchor <$> mws
+        memsn = joinHarmonium (average mzs) (mys >$< mxs) (average mws)
+     in (mprr,memsn,mtrns)
+
+-- | Direct expectation maximization for 'LatentProcess'es.
+latentProcessExpectationMaximization
+    :: ( ConjugatedLikelihood g x x w w, ConjugatedLikelihood f y x z w
+       , Transition Natural Mean w, Transition Natural Mean (AffineHarmonium g x x w w)
+       , Transition Mean Natural w
+       , Transition Mean Natural (AffineHarmonium f y x z w)
+       , Transition Mean Natural (AffineHarmonium g x x w w)
+       , Manifold (AffineHarmonium g x x w w)
+       , Bilinear g x x, Map Natural f x y, Bilinear f y x
+       , SamplePoint y ~ SamplePoint z )
+    => Observations (LatentProcess f g y x z w)
+    -> Natural # LatentProcess f g y x z w
+    -> Natural # LatentProcess f g y x z w
+latentProcessExpectationMaximization zss ltnt =
+    let (mprr,memsn,mtrns) = latentProcessExpectationStep zss ltnt
+        prr' = toNatural mprr
+        emsn' = fst . split $ toNatural memsn
+        trns' = fst . split $ toNatural mtrns
+     in joinLatentProcess prr' emsn' trns'
+
+-- | Expectation maximization for 'LatentProcess'es approximated through
+-- gradient ascent.
+latentProcessExpectationMaximizationAscent
+    :: ( ConjugatedLikelihood g x x w w, ConjugatedLikelihood f y x z w
+       , DuallyFlatExponentialFamily w
+       , LegendreExponentialFamily (AffineHarmonium f y x z w)
+       , LegendreExponentialFamily (AffineHarmonium g x x w w)
+       , Bilinear g x x, Map Natural f x y, Bilinear f y x
+       , SamplePoint y ~ SamplePoint z )
+    => Double
+    -> Int
+    -> GradientPursuit
+    -> [Sample z]
+    -> Natural # LatentProcess f g y x z w
+    -> Natural # LatentProcess f g y x z w
+latentProcessExpectationMaximizationAscent eps nstps gp zss ltnt =
+    let (mprr,mehrm,mthrm) = latentProcessExpectationStep zss ltnt
+        (nprr,nemsn,ntrns) = splitLatentProcess ltnt
+        neql0 = toNatural . snd $ split mehrm
+        neql1 = toNatural . snd $ split mthrm
+        nehrm = joinConjugatedHarmonium nemsn neql0
+        nthrm = joinConjugatedHarmonium ntrns neql1
+        nprr' = (!! nstps)
+            $ vanillaGradientSequence (relativeEntropyDifferential mprr) (-eps) gp nprr
+        nemsn' = fst . split . (!! nstps)
+            $ vanillaGradientSequence (relativeEntropyDifferential mehrm) (-eps) gp nehrm
+        ntrns' = fst . split . (!! nstps)
+            $ vanillaGradientSequence (relativeEntropyDifferential mthrm) (-eps) gp nthrm
+     in joinLatentProcess nprr' nemsn' ntrns'
diff --git a/Goal/Graphical/Models.hs b/Goal/Graphical/Models.hs
new file mode 100644
--- /dev/null
+++ b/Goal/Graphical/Models.hs
@@ -0,0 +1,43 @@
+{-# OPTIONS_GHC -fplugin=GHC.TypeLits.KnownNat.Solver -fplugin=GHC.TypeLits.Normalise -fconstraint-solver-iterations=10 #-}
+{-# LANGUAGE UndecidableInstances #-}
+
+-- | A few general definitions for Graphical Models.
+module Goal.Graphical.Models
+    ( Observation
+    , Observations
+    -- ** Hierarchical Models
+    , ObservablyContinuous ( logObservableDensities, observableDensities )
+    , logObservableDensity
+    , observableDensity
+    ) where
+
+--- Imports ---
+
+
+-- Package --
+
+import Goal.Geometry
+
+--- Latent Variable Class ---
+
+-- | An observation from a latent variable model.
+type family Observation f
+
+-- | A list of observations.
+type Observations f = [Observation f]
+
+-- | Probability densities over observations in a latent variable model.
+class ObservablyContinuous c f where
+    logObservableDensities :: c # f -> Observations f -> [Double]
+    logObservableDensities p = map (logObservableDensity p)
+
+    observableDensities :: c # f -> Observations f -> [Double]
+    observableDensities p = map (exp . logObservableDensity p)
+
+
+logObservableDensity :: ObservablyContinuous c f => c # f -> Observation f -> Double
+logObservableDensity p = head . logObservableDensities p . (:[])
+
+
+observableDensity :: ObservablyContinuous c f => c # f -> Observation f -> Double
+observableDensity p = exp . head . logObservableDensities p . (:[])
diff --git a/Goal/Graphical/Models/Dynamic.hs b/Goal/Graphical/Models/Dynamic.hs
new file mode 100644
--- /dev/null
+++ b/Goal/Graphical/Models/Dynamic.hs
@@ -0,0 +1,231 @@
+{-# OPTIONS_GHC -fplugin=GHC.TypeLits.KnownNat.Solver -fplugin=GHC.TypeLits.Normalise -fconstraint-solver-iterations=10 #-}
+{-# LANGUAGE UndecidableInstances #-}
+
+-- | 'Statistical' models where the observable biases depend on additional inputs.
+module Goal.Graphical.Models.Dynamic
+    (
+    LatentProcess (LatentProcess)
+    , HiddenMarkovModel
+    , SimpleKalmanFilter
+    , KalmanFilter
+    , sampleLatentProcess
+    -- ** Construction
+    , joinLatentProcess
+    , splitLatentProcess
+    -- ** Inference
+    , conjugatedFiltering
+    , conjugatedSmoothing
+    , conjugatedSmoothing0
+    ) where
+
+
+--- Imports  ---
+
+
+-- Goal --
+
+import Goal.Core
+import Goal.Geometry
+import Goal.Probability
+
+import Goal.Graphical.Models
+import Goal.Graphical.Inference
+import Goal.Graphical.Models.Harmonium
+
+import Data.List
+
+
+--- Generic ---
+
+
+-- | A conditional 'Harmonium', where the observable biases of the
+-- 'Harmonium' model depend on additional variables.
+newtype LatentProcess f g y x z w
+    = LatentProcess (AffineHarmonium f y x z w, Affine g x w x)
+
+type HiddenMarkovModel n k =
+    LatentProcess Tensor Tensor (Categorical n) (Categorical n) (Categorical k) (Categorical k)
+
+type SimpleKalmanFilter = LatentProcess Tensor Tensor NormalMean NormalMean Normal Normal
+
+type KalmanFilter n k
+  = LatentProcess Tensor Tensor (MVNMean n) (MVNMean k) (MultivariateNormal n) (MultivariateNormal k)
+
+type instance Observation (LatentProcess f g y x z w) = Sample z
+
+deriving instance (Manifold (AffineHarmonium f y x z w), Manifold (Affine g x w x))
+  => Manifold (LatentProcess f g y x z w)
+deriving instance (Manifold (AffineHarmonium f y x z w), Manifold (Affine g x w x))
+  => Product (LatentProcess f g y x z w)
+
+-- | Split a 'LatentProcess' into a prior, an emission distribution, and a
+-- transition distribution.
+splitLatentProcess
+    :: (Manifold z, Manifold w, Manifold (f y x), Manifold (g x x))
+    => c # LatentProcess f g y x z w
+    -> (c # w, c # Affine f y z x, c # Affine g x w x)
+splitLatentProcess ltnt =
+    let (hrm,trns) = split ltnt
+        (emsn,prr) = split hrm
+     in (prr,emsn,trns)
+
+-- | Construct a 'LatentProcess' from a prior, an emission distribution, and a
+-- transition distribution.
+joinLatentProcess
+    :: (Manifold z, Manifold w, Manifold (f y x), Manifold (g x x))
+    => c # w
+    -> c # Affine f y z x
+    -> c # Affine g x w x
+    -> c # LatentProcess f g y x z w
+joinLatentProcess prr emsn trns =
+    let hrm = join emsn prr
+     in join hrm trns
+
+latentProcessTransition
+    :: ( SamplePoint w ~ SamplePoint x, ExponentialFamily z
+       , Translation w x, Translation z y, Map Natural g x x
+       , ExponentialFamily x, Bilinear f x x
+       , Generative Natural w, Generative Natural z
+       , Bilinear g z x, Map Natural f y x )
+    => Natural # Affine f y z x -- ^ Emission Distribution
+    -> Natural # Affine g x w x -- ^ Transition Distribution
+    -> SamplePoint w
+    -> Random (SamplePoint (z,w))
+latentProcessTransition emsn trns w = do
+    w' <- samplePoint $ trns >.>* w
+    z' <- samplePoint $ emsn >.>* w'
+    return (z',w')
+
+-- | Generate a realization of the observable and latent states from a given
+-- latent process.
+sampleLatentProcess
+    :: ( SamplePoint w ~ SamplePoint x, ExponentialFamily z
+       , Translation w x, Translation z y, Map Natural g x x
+       , ExponentialFamily x, Bilinear f x x
+       , Generative Natural w, Generative Natural z
+       , Bilinear g z x, Map Natural f y x )
+    => Int
+    -> Natural # LatentProcess f g y x z w
+    -> Random (Sample (z,x))
+sampleLatentProcess n ltnt = do
+    let (prr,emsn,trns) = splitLatentProcess ltnt
+    x0 <- samplePoint prr
+    z0 <- samplePoint $ emsn >.>* x0
+    iterateM (n-1) (latentProcessTransition emsn trns . snd) (z0,x0)
+
+-- | Filtering for latent processes based on conjugated distributions.
+conjugatedFiltering
+    :: ( ConjugatedLikelihood g x x w w, Bilinear g x x
+       , ConjugatedLikelihood f y x z w, Bilinear f y x
+       , Map Natural g x x, Map Natural f x y )
+    => Natural # LatentProcess f g y x z w
+    -> Sample z
+    -> [Natural # w]
+conjugatedFiltering _ [] = []
+conjugatedFiltering ltnt (z:zs') =
+    let (prr,emsn,trns) = splitLatentProcess ltnt
+        prr' = conjugatedBayesRule emsn prr z
+     in scanl' (conjugatedForwardStep trns emsn) prr' zs'
+
+-- | Smoothing for latent processes based on conjugated distributions.
+conjugatedSmoothing
+    :: ( ConjugatedLikelihood g x x w w, Bilinear g x x
+       , ConjugatedLikelihood f y x z w, Bilinear f y x
+       , Map Natural g x x, Map Natural f x y )
+    => Natural # LatentProcess f g y x z w
+    -> Sample z
+    -> [Natural # w]
+conjugatedSmoothing ltnt zs =
+    let (prr,emsn,trns) = splitLatentProcess ltnt
+     in fst $ conjugatedSmoothing0 prr emsn trns zs
+
+-- | A more low-level implementation of smoothing which also returns joint
+-- distributions over current and subsequent states.
+conjugatedSmoothing0
+    :: ( ConjugatedLikelihood g x x w w, Bilinear g x x
+       , ConjugatedLikelihood f y x z w, Bilinear f y x
+       , Map Natural g x x, Map Natural f x y )
+    => Natural # w
+    -> Natural # Affine f y z x -- ^ Emission Distribution
+    -> Natural # Affine g x w x -- ^ Transition Distribution
+    -> Sample z
+    -> ([Natural # w],[Natural # AffineHarmonium g x x w w])
+conjugatedSmoothing0 _ _ _ [] = ([],[])
+conjugatedSmoothing0 prr emsn _ [z] =
+    ([conjugatedBayesRule emsn prr z],[])
+conjugatedSmoothing0 prr emsn trns (z:zs) =
+    let pst = conjugatedBayesRule emsn prr z
+        (trns',fwd) = splitConjugatedHarmonium . transposeHarmonium
+            $ joinConjugatedHarmonium trns pst
+        (smth:smths,hrms) = conjugatedSmoothing0 fwd emsn trns zs
+        hrm = transposeHarmonium $ joinConjugatedHarmonium trns' smth
+        bwd = snd $ splitConjugatedHarmonium hrm
+     in (bwd:smth:smths,hrm:hrms)
+
+
+--- Instances ---
+
+-- Implementations
+
+latentProcessLogDensity
+    :: ( ExponentialFamily z, ExponentialFamily x, Map Natural f y x
+       , Translation z y , Map Natural g x x, AbsolutelyContinuous Natural w
+       , SamplePoint w ~ SamplePoint x, AbsolutelyContinuous Natural z, Translation w x )
+    => Natural # w
+    -> Natural # Affine f y z x -- ^ Emission Distribution
+    -> Natural # Affine g x w x -- ^ Transition Distribution
+    -> Sample (z,w)
+    -> Double
+latentProcessLogDensity prr emsn trns zxs =
+    let (zs,xs) = unzip zxs
+        prrdns = logDensity prr $ head xs
+        trnsdnss = zipWith logDensity (trns >$>* xs) $ tail xs
+        emsndnss = zipWith logDensity (emsn >$>* xs) zs
+     in sum $ prrdns : trnsdnss ++ emsndnss
+
+conjugatedSmoothingLogDensity
+    :: ( ConjugatedLikelihood g x x w w, Bilinear g x x
+       , ConjugatedLikelihood f y x z w, Bilinear f y x
+       , Map Natural g x x, Map Natural f x y, ExponentialFamily y
+       , LegendreExponentialFamily z, LegendreExponentialFamily w )
+    => Natural # LatentProcess f g y x z w
+    -> Sample z
+    -> Double
+conjugatedSmoothingLogDensity ltnt zs =
+    let (_,emsn,_) = splitLatentProcess ltnt
+        smths = conjugatedSmoothing ltnt zs
+        hrms = joinConjugatedHarmonium emsn <$> smths
+     in sum $ zipWith logObservableDensity hrms zs
+
+-- Latent Processes
+
+instance Manifold (LatentProcess f g y x z w) => Statistical (LatentProcess f g y x z w) where
+    type SamplePoint (LatentProcess f g y x z w) = [SamplePoint (z,x)]
+
+instance ( ExponentialFamily z, ExponentialFamily x, Map Natural f y x
+         , Translation z y , Map Natural g x x, AbsolutelyContinuous Natural w
+         , SamplePoint w ~ SamplePoint x, AbsolutelyContinuous Natural z, Translation w x )
+  => AbsolutelyContinuous Natural (LatentProcess f g y x z w) where
+      logDensities ltnt zxss = do
+          zxs <- zxss
+          let (prr,emsn,trns) = splitLatentProcess ltnt
+          return $ latentProcessLogDensity prr emsn trns zxs
+
+instance ( ConjugatedLikelihood g x x w w, Bilinear g x x
+         , ConjugatedLikelihood f y x z w, Bilinear f y x
+         , Map Natural g x x, Map Natural f x y, ExponentialFamily y
+         , LegendreExponentialFamily z, LegendreExponentialFamily w )
+  => ObservablyContinuous Natural (LatentProcess f g y x z w) where
+    logObservableDensities ltnt = map (conjugatedSmoothingLogDensity ltnt)
+
+instance ( Manifold w , Manifold (g x x)
+         , Translation z y, Bilinear f y x )
+  => Translation (LatentProcess f g y x z w) y where
+    (>+>) ltnt y =
+        let (ehrm,trns) = split ltnt
+            (z,yx,w) = splitHarmonium ehrm
+            z' = z >+> y
+         in join (joinHarmonium z' yx w) trns
+    anchor ltnt =
+        anchor . snd . split . transposeHarmonium . fst $ split ltnt
+
diff --git a/Goal/Graphical/Models/Harmonium.hs b/Goal/Graphical/Models/Harmonium.hs
new file mode 100644
--- /dev/null
+++ b/Goal/Graphical/Models/Harmonium.hs
@@ -0,0 +1,687 @@
+{-# OPTIONS_GHC -fplugin=GHC.TypeLits.KnownNat.Solver -fplugin=GHC.TypeLits.Normalise -fconstraint-solver-iterations=10 #-}
+{-# LANGUAGE
+    TypeApplications,
+    UndecidableInstances,
+    NoStarIsType,
+    GeneralizedNewtypeDeriving,
+    StandaloneDeriving,
+    ScopedTypeVariables,
+    ExplicitNamespaces,
+    TypeOperators,
+    KindSignatures,
+    DataKinds,
+    RankNTypes,
+    TypeFamilies,
+    FlexibleContexts,
+    MultiParamTypeClasses,
+    ConstraintKinds,
+    FlexibleInstances
+#-}
+-- | An Exponential Family 'Harmonium' is a product exponential family with a
+-- particular bilinear structure (<https://papers.nips.cc/paper/2672-exponential-family-harmoniums-with-an-application-to-information-retrieval Welling, et al., 2005>).
+-- A 'Mixture' model is a special case of harmonium.
+module Goal.Graphical.Models.Harmonium
+    (
+    -- * Harmoniums
+      AffineHarmonium (AffineHarmonium)
+    , Harmonium
+    -- ** Constuction
+    , splitHarmonium
+    , joinHarmonium
+    -- ** Manipulation
+    , transposeHarmonium
+    -- ** Evaluation
+    , expectationStep
+    -- ** Sampling
+    , initialPass
+    , gibbsPass
+    -- ** Mixture Models
+    , Mixture
+    , AffineMixture
+    , joinNaturalMixture
+    , splitNaturalMixture
+    , joinMeanMixture
+    , splitMeanMixture
+    , joinSourceMixture
+    , splitSourceMixture
+    -- ** Linear Gaussian Harmoniums
+    , LinearGaussianHarmonium
+    -- ** Conjugated Harmoniums
+    , ConjugatedLikelihood (conjugationParameters)
+    , joinConjugatedHarmonium
+    , splitConjugatedHarmonium
+    ) where
+
+--- Imports ---
+
+
+import Goal.Core
+import Goal.Geometry
+import Goal.Probability
+
+import Goal.Graphical.Models
+
+import qualified Goal.Core.Vector.Storable as S
+
+
+--- Types ---
+
+
+-- | A 2-layer harmonium.
+newtype AffineHarmonium f y x z w = AffineHarmonium (Affine f y z x, w)
+
+deriving instance (Manifold z, Manifold (f y x), Manifold w)
+  => Manifold (AffineHarmonium f y x z w)
+deriving instance (Manifold z, Manifold (f y x), Manifold w)
+  => Product (AffineHarmonium f y x z w)
+
+type Harmonium f z w = AffineHarmonium f z w z w
+
+type instance Observation (AffineHarmonium f y x z w) = SamplePoint z
+
+-- | A 'Mixture' model is simply a 'AffineHarmonium' where the latent variable is
+-- 'Categorical'.
+type Mixture z k = Harmonium Tensor z (Categorical k)
+
+-- | A 'Mixture' where only a subset of the component parameters are mixed.
+type AffineMixture y z k =
+    AffineHarmonium Tensor y (Categorical k) z (Categorical k)
+
+type LinearGaussianHarmonium n k =
+    AffineHarmonium Tensor (MVNMean n) (MVNMean k) (MultivariateNormal n) (MultivariateNormal k)
+
+
+--- Classes ---
+
+
+-- | The conjugation parameters of a conjugated likelihood.
+class ( ExponentialFamily z, ExponentialFamily w, Map Natural f y x
+      , Translation z y , Translation w x
+      , SamplePoint y ~ SamplePoint z, SamplePoint x ~ SamplePoint w )
+  => ConjugatedLikelihood f y x z w where
+    conjugationParameters
+        :: Natural # Affine f y z x -- ^ Categorical likelihood
+        -> (Double, Natural # w) -- ^ Conjugation parameters
+
+
+--- Functions ---
+
+
+-- Construction --
+
+-- | Creates a 'Harmonium' from component parameters.
+joinHarmonium
+    :: (Manifold w, Manifold z, Manifold (f y x))
+    => c # z -- ^ Visible layer biases
+    -> c # f y x -- ^ ^ Interaction parameters
+    -> c # w -- ^ Hidden layer Biases
+    -> c # AffineHarmonium f y x z w -- ^ Harmonium
+joinHarmonium nz nyx = join (join nz nyx)
+
+-- | Splits a 'Harmonium' into component parameters.
+splitHarmonium
+    :: (Manifold z, Manifold (f y x), Manifold w)
+    => c # AffineHarmonium f y x z w -- ^ Harmonium
+    -> (c # z, c # f y x, c # w) -- ^ Biases and interaction parameters
+splitHarmonium hrm =
+    let (fzx,nw) = split hrm
+        (nz,nyx) = split fzx
+     in (nz,nyx,nw)
+
+-- | Build a mixture model in source coordinates.
+joinSourceMixture
+    :: (KnownNat k, Manifold z)
+    => S.Vector (k+1) (Source # z) -- ^ Mixture components
+    -> Source # Categorical k -- ^ Weights
+    -> Source # Mixture z k
+joinSourceMixture szs sx =
+    let (sz,szs') = S.splitAt szs
+        aff = join (S.head sz) (fromColumns szs')
+     in join aff sx
+
+-- | Build a mixture model in source coordinates.
+splitSourceMixture
+    :: (KnownNat k, Manifold z)
+    => Source # Mixture z k
+    -> (S.Vector (k+1) (Source # z), Source # Categorical k)
+splitSourceMixture mxmdl =
+    let (aff,sx) = split mxmdl
+        (sz0,szs0') = split aff
+     in (S.cons sz0 $ toColumns szs0' ,sx)
+
+-- | Build a mixture model in mean coordinates.
+joinMeanMixture
+    :: (KnownNat k, Manifold z)
+    => S.Vector (k+1) (Mean # z) -- ^ Mixture components
+    -> Mean # Categorical k -- ^ Weights
+    -> Mean # Mixture z k
+joinMeanMixture mzs mx =
+    let wghts = categoricalWeights mx
+        wmzs = S.zipWith (.>) wghts mzs
+        mz = S.foldr1 (+) wmzs
+        twmzs = S.tail wmzs
+        mzx = transpose . fromRows $ twmzs
+     in joinHarmonium mz mzx mx
+
+-- | Split a mixture model in mean coordinates.
+splitMeanMixture
+    :: ( KnownNat k, DuallyFlatExponentialFamily z )
+    => Mean # Mixture z k
+    -> (S.Vector (k+1) (Mean # z), Mean # Categorical k)
+splitMeanMixture hrm =
+    let (mz,mzx,mx) = splitHarmonium hrm
+        twmzs = toRows $ transpose mzx
+        wmzs = S.cons (mz - S.foldr (+) 0 twmzs) twmzs
+        wghts = categoricalWeights mx
+        mzs = S.zipWith (/>) wghts wmzs
+     in (mzs,mx)
+
+-- | A convenience function for building a categorical harmonium/mixture model.
+joinNaturalMixture
+    :: forall k z . ( KnownNat k, LegendreExponentialFamily z )
+    => S.Vector (k+1) (Natural # z) -- ^ Mixture components
+    -> Natural # Categorical k -- ^ Weights
+    -> Natural # Mixture z k -- ^ Mixture Model
+joinNaturalMixture nzs0 nx0 =
+    let nz0 :: S.Vector 1 (Natural # z)
+        (nz0,nzs0') = S.splitAt nzs0
+        nz = S.head nz0
+        nzs = S.map (subtract nz) nzs0'
+        nzx = fromMatrix . S.fromColumns $ S.map coordinates nzs
+        affzx = join nz nzx
+        rprms = snd $ conjugationParameters affzx
+        nx = nx0 - rprms
+     in joinHarmonium nz nzx nx
+
+-- | A convenience function for deconstructing a categorical harmonium/mixture model.
+splitNaturalMixture
+    :: forall k z . ( KnownNat k, LegendreExponentialFamily z )
+    => Natural # Mixture z k -- ^ Categorical harmonium
+    -> (S.Vector (k+1) (Natural # z), Natural # Categorical k) -- ^ (components, weights)
+splitNaturalMixture hrm =
+    let (nz,nzx,nx) = splitHarmonium hrm
+        affzx = join nz nzx
+        rprms = snd $ conjugationParameters affzx
+        nx0 = nx + rprms
+        nzs = S.map Point . S.toColumns $ toMatrix nzx
+        nzs0' = S.map (+ nz) nzs
+     in (S.cons nz nzs0',nx0)
+
+
+-- Manipulation --
+
+-- | Swap the biases and 'transpose' the interaction parameters of the given 'Harmonium'.
+transposeHarmonium
+    :: (Bilinear f y x, Manifold z, Manifold w)
+    => c # AffineHarmonium f y x z w
+    -> c # AffineHarmonium f x y w z
+transposeHarmonium hrm =
+        let (nz,nyx,nw) = splitHarmonium hrm
+         in joinHarmonium nw (transpose nyx) nz
+
+-- Evaluation --
+
+-- | Computes the joint expectations of a harmonium based on a sample from the
+-- observable layer.
+expectationStep
+    :: ( ExponentialFamily z, Map Natural f x y, Bilinear f y x
+       , Translation z y, Translation w x, LegendreExponentialFamily w )
+    => Sample z -- ^ Model Samples
+    -> Natural # AffineHarmonium f y x z w -- ^ Harmonium
+    -> Mean # AffineHarmonium f y x z w -- ^ Harmonium expected sufficient statistics
+expectationStep zs hrm =
+    let mzs = sufficientStatistic <$> zs
+        mys = anchor <$> mzs
+        pstr = fst . split $ transposeHarmonium hrm
+        mws = transition <$> pstr >$> mys
+        mxs = anchor <$> mws
+        myx = (>$<) mys mxs
+     in joinHarmonium (average mzs) myx $ average mws
+
+---- Sampling --
+
+-- | Initialize a Gibbs chain from a set of observations.
+initialPass
+    :: forall f x y z w
+    . ( ExponentialFamily z, Map Natural f x y, Manifold w
+      , SamplePoint y ~ SamplePoint z, Translation w x, Generative Natural w
+      , ExponentialFamily y, Bilinear f y x, LegendreExponentialFamily w )
+    => Natural # AffineHarmonium f y x z w -- ^ Harmonium
+    -> Sample z -- ^ Model Samples
+    -> Random (Sample (z, w))
+initialPass hrm zs = do
+    let pstr = fst . split $ transposeHarmonium hrm
+    ws <- mapM samplePoint $ pstr >$>* zs
+    return $ zip zs ws
+
+-- | Update a 'Sample' with Gibbs sampling.
+gibbsPass
+    :: ( ExponentialFamily z, Map Natural f x y, Translation z y
+       , Translation w x, SamplePoint z ~ SamplePoint y, Generative Natural w
+       , ExponentialFamily y, SamplePoint x ~ SamplePoint w, Bilinear f y x
+       , Map Natural f y x, ExponentialFamily x, Generative Natural z )
+    => Natural # AffineHarmonium f y x z w -- ^ Harmonium
+    -> Sample (z, w)
+    -> Random (Sample (z, w))
+gibbsPass hrm zws = do
+    let ws = snd <$> zws
+        pstr = fst . split $ transposeHarmonium hrm
+        lkl = fst $ split hrm
+    zs' <- mapM samplePoint $ lkl >$>* ws
+    ws' <- mapM samplePoint $ pstr >$>* zs'
+    return $ zip zs' ws'
+
+-- Conjugation --
+
+-- | The conjugation parameters of a conjugated `Harmonium`.
+harmoniumConjugationParameters
+    :: ConjugatedLikelihood f y x z w
+    => Natural # AffineHarmonium f y x z w -- ^ Categorical likelihood
+    -> (Double, Natural # w) -- ^ Conjugation parameters
+harmoniumConjugationParameters hrm =
+    conjugationParameters . fst $ split hrm
+
+-- | The conjugation parameters of a conjugated `Harmonium`.
+splitConjugatedHarmonium
+    :: ConjugatedLikelihood f y x z w
+    => Natural # AffineHarmonium f y x z w
+    -> (Natural # Affine f y z x, Natural # w) -- ^ Conjugation parameters
+splitConjugatedHarmonium hrm =
+    let (lkl,nw) = split hrm
+        cw = snd $ conjugationParameters lkl
+     in (lkl,nw + cw)
+
+-- | The conjugation parameters of a conjugated `Harmonium`.
+joinConjugatedHarmonium
+    :: ConjugatedLikelihood f y x z w
+    => Natural # Affine f y z x -- ^ Conjugation parameters
+    -> Natural # w
+    -> Natural # AffineHarmonium f y x z w -- ^ Categorical likelihood
+joinConjugatedHarmonium lkl nw =
+    let cw = snd $ conjugationParameters lkl
+     in join lkl $ nw - cw
+
+-- | The conjugation parameters of a conjugated `Harmonium`.
+sampleConjugated
+    :: forall f y x z w
+     . ( ConjugatedLikelihood f y x z w, Generative Natural w
+       , Generative Natural z, Map Natural f y x )
+    => Int
+    -> Natural # AffineHarmonium f y x z w -- ^ Categorical likelihood
+    -> Random (Sample (z,w)) -- ^ Conjugation parameters
+sampleConjugated n hrm = do
+    let (lkl,nw) = split hrm
+        nw' = nw + snd (conjugationParameters lkl)
+    ws <- sample n nw'
+    let mws :: [Mean # w]
+        mws = sufficientStatistic <$> ws
+    zs <- mapM samplePoint $ lkl >$+> mws
+    return $ zip zs ws
+
+-- | The conjugation parameters of a conjugated `Harmonium`.
+conjugatedPotential
+    :: ( LegendreExponentialFamily w, ConjugatedLikelihood f y x z w )
+    => Natural # AffineHarmonium f y x z w -- ^ Categorical likelihood
+    -> Double -- ^ Conjugation parameters
+conjugatedPotential hrm = do
+    let (lkl,nw) = split hrm
+        (rho0,rprms) = conjugationParameters lkl
+     in potential (nw + rprms) + rho0
+
+
+--- Internal ---
+
+
+-- Conjugation --
+
+-- | The unnormalized density of a given 'Harmonium' 'Point'.
+unnormalizedHarmoniumObservableLogDensity
+    :: forall f y x z w
+    . ( ExponentialFamily z, ExponentialFamily y
+      , LegendreExponentialFamily w, Translation w x, Translation z y
+      , Map Natural f x y, Bilinear f y x )
+    => Natural # AffineHarmonium f y x z w
+    -> Sample z
+    -> [Double]
+unnormalizedHarmoniumObservableLogDensity hrm zs =
+    let (pstr,nz) = split $ transposeHarmonium hrm
+        mzs = sufficientStatistic <$> zs
+        nrgs = zipWith (+) (dotMap nz mzs) $ potential <$> pstr >$+> mzs
+     in zipWith (+) nrgs $ logBaseMeasure (Proxy @ z) <$> zs
+
+--- | Computes the negative log-likelihood of a sample point of a conjugated harmonium.
+logConjugatedDensities
+    :: forall f x y z w
+    . ( Bilinear f y x, Translation z y
+      , LegendreExponentialFamily z, ExponentialFamily y
+      , LegendreExponentialFamily w, Translation w x, Map Natural f x y)
+      => (Double, Natural # w) -- ^ Conjugation Parameters
+      -> Natural # AffineHarmonium f y x z w
+      -> Sample z
+      -> [Double]
+logConjugatedDensities (rho0,rprms) hrm z =
+    let udns = unnormalizedHarmoniumObservableLogDensity hrm z
+        nx = snd $ split hrm
+     in subtract (potential (nx + rprms) + rho0) <$> udns
+
+-- Mixtures --
+
+mixtureLikelihoodConjugationParameters
+    :: (KnownNat k, LegendreExponentialFamily z, Translation z y)
+    => Natural # Affine Tensor y z (Categorical k) -- ^ Categorical likelihood
+    -> (Double, Natural # Categorical k) -- ^ Conjugation parameters
+mixtureLikelihoodConjugationParameters aff =
+    let (nz,nyx) = split aff
+        rho0 = potential nz
+        rprms = S.map (\nyxi -> subtract rho0 . potential $ nz >+> nyxi) $ toColumns nyx
+     in (rho0, Point rprms)
+
+affineMixtureToMixture
+    :: (KnownNat k, Manifold z, Manifold y, Translation z y)
+    => Natural # AffineMixture y z k
+    -> Natural # Mixture z k
+affineMixtureToMixture lmxmdl =
+    let (flsk,nk) = split lmxmdl
+        (nls,nlk) = split flsk
+        nlsk = fromColumns . S.map (0 >+>) $ toColumns nlk
+     in join (join nls nlsk) nk
+
+mixtureToAffineMixture
+    :: (KnownNat k, Manifold y, Manifold z, Translation z y)
+    => Mean # Mixture z k
+    -> Mean # AffineMixture y z k
+mixtureToAffineMixture mxmdl =
+    let (flsk,mk) = split mxmdl
+        (mls,mlsk) = split flsk
+        mlk = fromColumns . S.map anchor $ toColumns mlsk
+     in join (join mls mlk) mk
+
+
+-- Linear Gaussian Harmoniums --
+
+linearGaussianHarmoniumConjugationParameters
+    :: (KnownNat n, KnownNat k)
+    => Natural # Affine Tensor (MVNMean n) (MultivariateNormal n) (MVNMean k)
+    -> (Double, Natural # MultivariateNormal k) -- ^ Conjugation parameters
+linearGaussianHarmoniumConjugationParameters aff =
+    let (thts,tht30) = split aff
+        (tht1,tht2) = splitNaturalMultivariateNormal thts
+        tht3 = toMatrix tht30
+        ttht3 = S.transpose tht3
+        itht2 = S.pseudoInverse tht2
+        rho0 = -0.25 * tht1 `S.dotProduct` (itht2 `S.matrixVectorMultiply` tht1)
+            -0.5 * (log . S.determinant . negate $ 2*tht2)
+        rho1 = -0.5 * ttht3 `S.matrixVectorMultiply` (itht2 `S.matrixVectorMultiply` tht1)
+        rho2 = -0.25 * ttht3 `S.matrixMatrixMultiply` (itht2 `S.matrixMatrixMultiply` tht3)
+     in (rho0, joinNaturalMultivariateNormal rho1 rho2)
+
+univariateToLinearGaussianHarmonium
+    :: c # AffineHarmonium Tensor NormalMean NormalMean Normal Normal
+    -> c # LinearGaussianHarmonium 1 1
+univariateToLinearGaussianHarmonium hrm =
+    let (z,zx,x) = splitHarmonium hrm
+     in joinHarmonium (breakPoint z) (breakPoint zx) (breakPoint x)
+
+linearGaussianHarmoniumToUnivariate
+    :: c # LinearGaussianHarmonium 1 1
+    -> c # AffineHarmonium Tensor NormalMean NormalMean Normal Normal
+linearGaussianHarmoniumToUnivariate hrm =
+    let (z,zx,x) = splitHarmonium hrm
+     in joinHarmonium (breakPoint z) (breakPoint zx) (breakPoint x)
+
+univariateToLinearModel
+    :: Natural # Affine Tensor NormalMean Normal NormalMean
+    -> Natural # Affine Tensor (MVNMean 1) (MultivariateNormal 1) (MVNMean 1)
+univariateToLinearModel aff =
+    let (z,zx) = split aff
+     in join (breakPoint z) (breakPoint zx)
+
+naturalLinearGaussianHarmoniumToJoint
+    :: (KnownNat n, KnownNat k)
+    => Natural # LinearGaussianHarmonium n k
+    -> Natural # MultivariateNormal (n+k)
+naturalLinearGaussianHarmoniumToJoint hrm =
+    let (z,zx,x) = splitHarmonium hrm
+        zxmtx = toMatrix zx/2
+        mvnz = splitNaturalMultivariateNormal z
+        mvnx = splitNaturalMultivariateNormal x
+        (mu,cvr) = fromLinearGaussianHarmonium0 mvnz zxmtx mvnx
+     in joinNaturalMultivariateNormal mu cvr
+
+naturalJointToLinearGaussianHarmonium
+    :: (KnownNat n, KnownNat k)
+    => Natural # MultivariateNormal (n+k)
+    -> Natural # LinearGaussianHarmonium n k
+naturalJointToLinearGaussianHarmonium mvn =
+    let (mu,cvr) = splitNaturalMultivariateNormal mvn
+        ((muz,cvrz),zxmtx,(mux,cvrx)) = toLinearGaussianHarmonium0 mu cvr
+        zx = 2*fromMatrix zxmtx
+        z = joinNaturalMultivariateNormal muz cvrz
+        x = joinNaturalMultivariateNormal mux cvrx
+     in joinHarmonium z zx x
+
+meanLinearGaussianHarmoniumToJoint
+    :: (KnownNat n, KnownNat k)
+    => Mean # LinearGaussianHarmonium n k
+    -> Mean # MultivariateNormal (n+k)
+meanLinearGaussianHarmoniumToJoint hrm =
+    let (z,zx,x) = splitHarmonium hrm
+        zxmtx = toMatrix zx
+        mvnz = splitMeanMultivariateNormal z
+        mvnx = splitMeanMultivariateNormal x
+        (mu,cvr) = fromLinearGaussianHarmonium0 mvnz zxmtx mvnx
+     in joinMeanMultivariateNormal mu cvr
+
+meanJointToLinearGaussianHarmonium
+    :: (KnownNat n, KnownNat k)
+    => Mean # MultivariateNormal (n+k)
+    -> Mean # LinearGaussianHarmonium n k
+meanJointToLinearGaussianHarmonium mvn =
+    let (mu,cvr) = splitMeanMultivariateNormal mvn
+        ((muz,cvrz),zxmtx,(mux,cvrx)) = toLinearGaussianHarmonium0 mu cvr
+        zx = fromMatrix zxmtx
+        z = joinMeanMultivariateNormal muz cvrz
+        x = joinMeanMultivariateNormal mux cvrx
+     in joinHarmonium z zx x
+
+fromLinearGaussianHarmonium0
+    :: (KnownNat n, KnownNat k)
+    => (S.Vector n Double, S.Matrix n n Double)
+    -> S.Matrix n k Double
+    -> (S.Vector k Double, S.Matrix k k Double)
+    -> (S.Vector (n+k) Double, S.Matrix (n+k) (n+k) Double)
+fromLinearGaussianHarmonium0 (muz,cvrz) zxmtx (mux,cvrx) =
+    let mu = muz S.++ mux
+        top = S.horizontalConcat cvrz zxmtx
+        btm = S.horizontalConcat (S.transpose zxmtx) cvrx
+     in (mu, S.verticalConcat top btm)
+
+toLinearGaussianHarmonium0
+    :: (KnownNat n, KnownNat k)
+    => S.Vector (n+k) Double
+    -> S.Matrix (n+k) (n+k) Double
+    -> ( (S.Vector n Double, S.Matrix n n Double)
+       , S.Matrix n k Double
+       , (S.Vector k Double, S.Matrix k k Double) )
+toLinearGaussianHarmonium0 mu cvr =
+    let (muz,mux) = S.splitAt mu
+        (tops,btms) = S.splitAt $ S.toRows cvr
+        (cvrzs,zxmtxs) = S.splitAt . S.toColumns $ S.fromRows tops
+        cvrz = S.fromColumns cvrzs
+        zxmtx = S.fromColumns zxmtxs
+        cvrx = S.fromColumns . S.drop . S.toColumns $ S.fromRows btms
+     in ((muz,cvrz),zxmtx,(mux,cvrx))
+
+harmoniumLogBaseMeasure
+    :: forall f y x z w . (ExponentialFamily z, ExponentialFamily w)
+    => Proxy (AffineHarmonium f y x z w)
+    -> SamplePoint (z,w)
+    -> Double
+harmoniumLogBaseMeasure _ (z,w) =
+    logBaseMeasure (Proxy @ z) z + logBaseMeasure (Proxy @ w) w
+
+
+--- Instances ---
+
+
+instance Manifold (AffineHarmonium f y x z w) => Statistical (AffineHarmonium f y x z w) where
+    type SamplePoint (AffineHarmonium f y x z w) = SamplePoint (z,w)
+
+instance ( ExponentialFamily z, ExponentialFamily x, Translation z y
+         , Translation w x
+         , ExponentialFamily w, ExponentialFamily y, Bilinear f y x )
+  => ExponentialFamily (AffineHarmonium f y x z w) where
+      sufficientStatistic (z,w) =
+          let mz = sufficientStatistic z
+              mw = sufficientStatistic w
+              my = anchor mz
+              mx = anchor mw
+           in joinHarmonium mz (my >.< mx) mw
+      averageSufficientStatistic zws =
+          let (zs,ws) = unzip zws
+              mzs = sufficientStatistic <$> zs
+              mws = sufficientStatistic <$> ws
+              mys = anchor <$> mzs
+              mxs = anchor <$> mws
+           in joinHarmonium (average mzs) (mys >$< mxs) (average mws)
+      logBaseMeasure = harmoniumLogBaseMeasure
+
+instance ( KnownNat k, LegendreExponentialFamily z
+         , Translation y z, LegendreExponentialFamily y, SamplePoint z ~ SamplePoint y)
+  => ConjugatedLikelihood Tensor z (Categorical k) y (Categorical k) where
+    conjugationParameters = mixtureLikelihoodConjugationParameters
+
+instance ConjugatedLikelihood Tensor NormalMean NormalMean Normal Normal where
+    conjugationParameters aff =
+        let rprms :: Natural # MultivariateNormal 1
+            (rho0,rprms) = conjugationParameters $ univariateToLinearModel aff
+         in (rho0,breakPoint rprms)
+
+instance (KnownNat n, KnownNat k) => ConjugatedLikelihood Tensor (MVNMean n) (MVNMean k)
+    (MultivariateNormal n) (MultivariateNormal k) where
+        conjugationParameters = linearGaussianHarmoniumConjugationParameters
+
+--instance ( KnownNat k, LegendreExponentialFamily z
+--         , Generative Natural z, Manifold (Mixture z k) )
+--         => Generative Natural (Mixture z k) where
+--    sample = sampleConjugated
+
+--instance (KnownNat k, LegendreExponentialFamily z)
+--  => Transition Natural Mean (Mixture z k) where
+--    transition nhrm =
+--        let (nzs,nx) = splitNaturalMixture nhrm
+--            mx = toMean nx
+--            mzs = S.map transition nzs
+--         in joinMeanMixture mzs mx
+
+instance ( KnownNat k, Manifold y, Manifold z
+         , LegendreExponentialFamily z, Translation z y )
+  => Transition Natural Mean (AffineMixture y z k) where
+    transition mxmdl0 =
+        let mxmdl = affineMixtureToMixture mxmdl0
+            (nzs,nx) = splitNaturalMixture mxmdl
+            mx = toMean nx
+            mzs = S.map transition nzs
+         in mixtureToAffineMixture $ joinMeanMixture mzs mx
+
+instance ( KnownNat k, Manifold y, Manifold z, LegendreExponentialFamily z
+         , Generative Natural z, Translation z y )
+  => Generative Natural (AffineMixture y z k) where
+      sample n = sampleConjugated n . affineMixtureToMixture
+
+instance (KnownNat k, DuallyFlatExponentialFamily z)
+  => Transition Mean Natural (Mixture z k) where
+    transition mhrm =
+        let (mzs,mx) = splitMeanMixture mhrm
+            nx = transition mx
+            nzs = S.map transition mzs
+         in joinNaturalMixture nzs nx
+
+instance (KnownNat k, LegendreExponentialFamily z, Transition Natural Source z)
+  => Transition Natural Source (Mixture z k) where
+    transition nhrm =
+        let (nzs,nx) = splitNaturalMixture nhrm
+            sx = transition nx
+            szs = S.map transition nzs
+         in joinSourceMixture szs sx
+
+instance (KnownNat k, LegendreExponentialFamily z, Transition Source Natural z)
+  => Transition Source Natural (Mixture z k) where
+    transition shrm =
+        let (szs,sx) = splitSourceMixture shrm
+            nx = transition sx
+            nzs = S.map transition szs
+         in joinNaturalMixture nzs nx
+
+instance Transition Natural Mean
+  (AffineHarmonium Tensor NormalMean NormalMean Normal Normal) where
+      transition = linearGaussianHarmoniumToUnivariate . transition . univariateToLinearGaussianHarmonium
+
+instance Transition Mean Natural
+  (AffineHarmonium Tensor NormalMean NormalMean Normal Normal) where
+      transition =  linearGaussianHarmoniumToUnivariate . transition . univariateToLinearGaussianHarmonium
+
+instance (KnownNat n, KnownNat k) => Transition Natural Mean
+  (AffineHarmonium Tensor (MVNMean n) (MVNMean k)
+    (MultivariateNormal n) (MultivariateNormal k)) where
+      transition = meanJointToLinearGaussianHarmonium . transition
+        . naturalLinearGaussianHarmoniumToJoint
+
+instance (KnownNat n, KnownNat k) => Transition Mean Natural
+  (AffineHarmonium Tensor (MVNMean n) (MVNMean k)
+    (MultivariateNormal n) (MultivariateNormal k)) where
+      transition = naturalJointToLinearGaussianHarmonium . transition
+        . meanLinearGaussianHarmoniumToJoint
+
+--type instance PotentialCoordinates (Mixture z k) = Natural
+--
+--instance (KnownNat k, LegendreExponentialFamily z) => Legendre (Mixture z k) where
+--      potential = conjugatedPotential
+
+type instance PotentialCoordinates (AffineHarmonium f y x z w) = Natural
+
+instance ( Manifold (f y x), LegendreExponentialFamily w, ConjugatedLikelihood f y x z w )
+  => Legendre (AffineHarmonium f y x z w) where
+      potential = conjugatedPotential
+
+instance ( Manifold (f y x), LegendreExponentialFamily w
+         , Transition Mean Natural (AffineHarmonium f y x z w), ConjugatedLikelihood f y x z w )
+  => DuallyFlat (AffineHarmonium f y x z w) where
+    dualPotential mhrm =
+        let nhrm = toNatural mhrm
+         in mhrm <.> nhrm - potential nhrm
+
+instance ( Bilinear f y x, ExponentialFamily y, ExponentialFamily x
+         , LegendreExponentialFamily w, ConjugatedLikelihood f y x z w )
+  => AbsolutelyContinuous Natural (AffineHarmonium f y x z w) where
+    logDensities = exponentialFamilyLogDensities
+
+instance ( ConjugatedLikelihood f y x z w, LegendreExponentialFamily z
+         , ExponentialFamily y, LegendreExponentialFamily w
+         , Map Natural f x y, Bilinear f x y )
+  => ObservablyContinuous Natural (AffineHarmonium f y x z w) where
+    logObservableDensities hrm zs =
+        let rho0rprms = harmoniumConjugationParameters hrm
+         in logConjugatedDensities rho0rprms hrm zs
+
+instance ( LegendreExponentialFamily z, LegendreExponentialFamily w
+         , ConjugatedLikelihood f y x z w, Map Natural f x y, Bilinear f x y
+         , LegendreExponentialFamily (AffineHarmonium f y x z w)
+         , Manifold (f y x), SamplePoint z ~ t, ExponentialFamily y)
+  => LogLikelihood Natural (AffineHarmonium f y x z w) t where
+    logLikelihood xs hrm =
+         average $ logObservableDensities hrm xs
+    logLikelihoodDifferential zs hrm =
+        let pxs = expectationStep zs hrm
+            qxs = transition hrm
+         in pxs - qxs
+
+instance ( Translation z y, Manifold w, Manifold (f y x) )
+  => Translation (AffineHarmonium f y x z w) y where
+      (>+>) hrm ny =
+          let (nz,nyx,nw) = splitHarmonium hrm
+           in joinHarmonium (nz >+> ny) nyx nw
+      anchor hrm =
+          let (nz,_,_) = splitHarmonium hrm
+           in anchor nz
+
diff --git a/Goal/Graphical/Models/Harmonium/FactorAnalysis.hs b/Goal/Graphical/Models/Harmonium/FactorAnalysis.hs
new file mode 100644
--- /dev/null
+++ b/Goal/Graphical/Models/Harmonium/FactorAnalysis.hs
@@ -0,0 +1,107 @@
+{-# OPTIONS_GHC -fplugin=GHC.TypeLits.KnownNat.Solver -fplugin=GHC.TypeLits.Normalise -fconstraint-solver-iterations=10 #-}
+{-# LANGUAGE
+    TypeApplications,
+    UndecidableInstances,
+    NoStarIsType,
+    GeneralizedNewtypeDeriving,
+    StandaloneDeriving,
+    ScopedTypeVariables,
+    ExplicitNamespaces,
+    TypeOperators,
+    KindSignatures,
+    DataKinds,
+    RankNTypes,
+    TypeFamilies,
+    FlexibleContexts,
+    MultiParamTypeClasses,
+    ConstraintKinds,
+    FlexibleInstances
+#-}
+-- | An Exponential Family 'Harmonium' is a product exponential family with a
+-- particular bilinear structure (<https://papers.nips.cc/paper/2672-exponential-family-harmoniums-with-an-application-to-information-retrieval Welling, et al., 2005>).
+-- A 'Mixture' model is a special case of harmonium. A 'FactorAnalysis' model
+-- can also be interpreted as a 'Harmonium' with a fixed latent distribution.
+module Goal.Graphical.Models.Harmonium.FactorAnalysis
+    (
+    -- * Factor Analysis
+      FactorAnalysis
+    , factorAnalysisObservableDistribution
+    , factorAnalysisExpectationMaximization
+    , factorAnalysisUniqueness
+    ) where
+
+--- Imports ---
+
+
+import Goal.Core
+import Goal.Geometry
+import Goal.Probability
+
+import Goal.Graphical.Models
+import Goal.Graphical.Models.Harmonium
+
+import qualified Goal.Core.Vector.Storable as S
+
+
+--- Types ---
+
+
+type FactorAnalysis n k = Affine Tensor (MVNMean n) (Replicated n Normal) (MVNMean k)
+
+type instance Observation (FactorAnalysis n k) = S.Vector n Double
+
+factorAnalysisObservableDistribution
+    :: (KnownNat n, KnownNat k)
+    => Natural # FactorAnalysis n k
+    -> Natural # MultivariateNormal n
+factorAnalysisObservableDistribution =
+     snd . splitConjugatedHarmonium . transposeHarmonium
+     . naturalFactorAnalysisToLGH
+
+factorAnalysisExpectationMaximization
+    :: ( KnownNat n, KnownNat k)
+    => [S.Vector n Double]
+    -> Natural # FactorAnalysis n k
+    -> Natural # FactorAnalysis n k
+factorAnalysisExpectationMaximization zs fa =
+    transition .sourceFactorAnalysisMaximizationStep . expectationStep zs
+        $ naturalFactorAnalysisToLGH fa
+
+factorAnalysisUniqueness
+    :: (KnownNat n, KnownNat k)
+    => Natural # FactorAnalysis n k
+    -> S.Vector n Double
+factorAnalysisUniqueness fa =
+    let lds = toMatrix . snd . split $ toSource fa
+        sgs = S.takeDiagonal . snd . splitMultivariateNormal . toSource $ factorAnalysisObservableDistribution fa
+        cms = S.takeDiagonal . S.matrixMatrixMultiply lds $ S.transpose lds
+     in (sgs - cms) / sgs
+
+
+-- Internal --
+
+naturalFactorAnalysisToLGH
+    :: (KnownNat n, KnownNat k)
+    => Natural # FactorAnalysis n k
+    -> Natural # LinearGaussianHarmonium n k
+naturalFactorAnalysisToLGH fa =
+    let (nzs,tns) = split fa
+        (mus,vrs) = S.toPair . S.toColumns . S.fromRows . S.map coordinates $ splitReplicated nzs
+        mvn = joinNaturalMultivariateNormal mus $ S.diagonalMatrix vrs
+        fa' = join mvn tns
+     in joinConjugatedHarmonium fa' $ toNatural . joinMultivariateNormal 0 $ S.diagonalMatrix 1
+
+sourceFactorAnalysisMaximizationStep
+    :: forall n k . (KnownNat n, KnownNat k)
+    => Mean # LinearGaussianHarmonium n k
+    -> Source # FactorAnalysis n k
+sourceFactorAnalysisMaximizationStep hrm =
+    let (nz,nzx,nx) = splitHarmonium hrm
+        (muz,etaz) = splitMeanMultivariateNormal nz
+        (mux,etax) = splitMeanMultivariateNormal nx
+        outrs = toMatrix nzx - S.outerProduct muz mux
+        wmtx = S.matrixMatrixMultiply outrs $ S.inverse etax
+        zcvr = etaz - S.outerProduct muz muz
+        vrs = S.takeDiagonal $ zcvr - S.matrixMatrixMultiply wmtx (S.transpose outrs)
+        snrms = joinReplicated $ S.zipWith (curry fromTuple) muz vrs
+     in join snrms $ fromMatrix wmtx
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/README.md b/README.md
new file mode 100644
--- /dev/null
+++ b/README.md
@@ -0,0 +1,66 @@
+This library provides definitions and algorithms for various graphical models
+such as mixture models, kalman Filters, and restricted Boltzmann machines, as
+well as algorithms for fitting them e.g.  expectation maximization and
+contrastive divergence minimization. Underlying all of these models is is a
+generalized linear object known as a Harmonium, and in the following I will
+briefly introduce them.
+
+The core definition of this library is a `Manifold` of joint distributions I
+call of an `AffineHarmonium`
+```haskell
+newtype AffineHarmonium f y x z w = AffineHarmonium (Affine f y z x, w)
+```
+which is a product `Manifold` composed of a of a `Manifold` of likelihood
+functions `Affine f y z x`, and a `Manifold` of distributions `w` that partially
+define the space of priors. `AffineHarmoniums` provide a bit more flexibility
+than what I call a `Harmonium`
+```haskell
+type Harmonium f z w = AffineHarmonium f z w z w
+```
+which is a simpler object. Nevertheless, from a theoretical point of view, an
+`AffineHarmonium` is a special case of a `Harmonium`, and we may think of them
+more or less equivalently.
+
+A `Harmonium` is a model over a observable variables and latent variables, and represents a sort of generalized linear joint distribution over the two of them. The theory for `Harmonium`s is summarized well by [this paper](https://papers.nips.cc/paper/2004/hash/0e900ad84f63618452210ab8baae0218-Abstract.html) and [this paper](https://proceedings.neurips.cc/paper/2013/hash/28f0b864598a1291557bed248a998d4e-Abstract.html). Although `Harmonium`s might seem like a little-studied, and esoteric object, various well known models, such as mixture models and restricted Boltzmann machines, are in fact `Harmonium`s, and various other models, such as factor analysis, Kalman filters, or hidden Markov models, can be expressed in terms of them.
+
+All of the aforementioned models can be fit with Expectation-Maximization (EM), and EM can be expressed in an entirely general manner for `Harmonium`s. Firstly, the expectation step of a `Harmonium` is implemented by
+```haskell
+expectationStep
+    :: ( ExponentialFamily z, Map Natural f x y, Bilinear f y x
+       , Translation z y, Translation w x, LegendreExponentialFamily w )
+    => Sample z -- ^ Model Samples
+    -> Natural # AffineHarmonium f y x z w -- ^ Harmonium
+    -> Mean # AffineHarmonium f y x z w -- ^ Harmonium expected sufficient statistics
+expectationStep zs hrm =
+    let mzs = sufficientStatistic <$> zs
+        mys = anchor <$> mzs
+        pstr = fst . split $ transposeHarmonium hrm
+        mws = transition <$> pstr >$> mys
+        mxs = anchor <$> mws
+        myx = (>$<) mys mxs
+     in joinHarmonium (average mzs) myx $ average mws
+```
+In summary, what we do is
+
+- take some observations,
+- compute their `sufficientStatistics`,
+- map these statistics into the predictions of the latent variables,
+- transition these latent predictions from `Natural` coordinates to `Mean` coordinates,
+- and assemble the results into the `Mean` `sufficientStatistics` of the joint distribution.
+
+The maximization step then consists simply of mapping the whole joint distribution from `Mean` back to `Natural` coordinates, such that `expectationMaximization` may be expressed as
+```haskell
+expectationMaximization
+    :: ( DuallyFlatExponentialFamily (AffineHarmonium f y x z w)
+       , ExponentialFamily z, Map Natural f x y, Bilinear f y x
+       , Translation z y, Translation w x, LegendreExponentialFamily w )
+    => Sample z
+    -> Natural # AffineHarmonium f y x z w
+    -> Natural # AffineHarmonium f y x z w
+expectationMaximization zs hrm = transition $ expectationStep zs hrm
+```
+
+As such, for a wide variety of models, we may reduce implementing expectation maximization to instantating the class requirements of the `expectationMaximization` function. This is rarely trivial, but in some sense, much more straight-forward and well-defined that deriving EM algorithms from scratch.
+
+For in-depth tutorials visit my
+[blog](https://sacha-sokoloski.gitlab.io/website/pages/blog.html).
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/benchmarks/com-poisson.hs b/benchmarks/com-poisson.hs
new file mode 100644
--- /dev/null
+++ b/benchmarks/com-poisson.hs
@@ -0,0 +1,119 @@
+{-# LANGUAGE DataKinds,ScopedTypeVariables,FlexibleContexts,TypeOperators,TypeFamilies #-}
+
+import Goal.Geometry
+import Goal.Probability
+import Goal.Graphical
+
+import qualified Criterion.Main as C
+
+
+--- Globals ---
+
+
+type N = 101
+type K = 10
+
+type PoissonMixture = Mixture (Replicated N Poisson) K
+type CoMPoissonMixture = AffineMixture (Replicated N Poisson) (Replicated N CoMPoisson) K
+
+mixtureStep
+    :: ( LegendreExponentialFamily w, ConjugatedLikelihood f y x z w
+       , ExponentialFamily x, ExponentialFamily y
+       , Transition Natural Mean (AffineHarmonium f y x z w)
+       , Bilinear f y x, Map Natural f x y )
+       => [SamplePoint z]
+       -> (Natural # AffineHarmonium f y x z w)
+       -> Natural # AffineHarmonium f y x z w
+mixtureStep zs mxmdl =
+    expectationMaximizationAscent 2e-3 defaultAdamPursuit zs mxmdl !! 20
+
+comDeviations
+    :: Double
+    -> Double
+    -> Mean # CoMPoisson
+    -> Natural # CoMPoisson
+    -> (Double,Double,Double)
+comDeviations err lgprt0 mcm0 ncm =
+    let lgprt = comPoissonLogPartitionSum err ncm
+        mcm = comPoissonMeans err ncm
+        [muerr,nuerr] = listCoordinates $ mcm0 - mcm
+     in (lgprt0 - lgprt, muerr,nuerr)
+
+comSDs
+    :: [Double]
+    -> [Mean # CoMPoisson]
+    -> [Natural # CoMPoisson]
+    -> Double
+    -> (Double,Double,Double)
+comSDs lgprts0 mcm0s ncms err =
+    let (lgprterrs,muerrs,nuerrs) = unzip3 $ zipWith3 (comDeviations err) lgprts0 mcm0s ncms
+        lgprtsd = snd $ estimateMeanVariance lgprterrs
+        musd = snd $ estimateMeanVariance muerrs
+        nusd = snd $ estimateMeanVariance nuerrs
+     in (lgprtsd,musd,nusd)
+
+comPoissonMeans :: Double -> Natural # CoMPoisson -> Mean # CoMPoisson
+comPoissonMeans eps cp =
+    let ss :: Int -> Mean # CoMPoisson
+        ss = sufficientStatistic
+     in Point $ comPoissonExpectations eps (coordinates . ss) cp
+
+
+
+--- Main ---
+
+
+main :: IO ()
+main = do
+
+    mx0 :: Natural # PoissonMixture
+        <- realize $ uniformInitialize (-2,2)
+
+    zxs <- realize $ sample 100 mx0
+    let zs = fst <$> zxs
+
+    let (nyx,ny) = split $ transposeHarmonium mx0
+        cmx0 :: Natural # CoMPoissonMixture
+        cmx0 = transposeHarmonium . join nyx
+            $ mapReplicatedPoint (`join` (-1.5)) ny
+
+    let cms :: [Source # CoMPoisson]
+        cms = do
+            mu <- [0.5,2,20]
+            nu <- [0.3,1,10]
+            return $ fromTuple (mu,nu)
+        ncms = toNatural <$> cms
+
+    let err0 = 1e-20
+    let lgprt0s = comPoissonLogPartitionSum err0 <$> ncms
+        mcm0s = comPoissonMeans err0 <$> ncms
+
+    let errs = [1e-2,1e-4,1e-6,1e-8,1e-10,1e-12,1e-14,1e-16]
+        comFun = comSDs lgprt0s mcm0s ncms
+
+    sequence_ $ do
+        err <- errs
+        return $ do
+            putStrLn $ concat ["Error: ", show err]
+            putStrLn $ concat ["Deviations: ", show $ comFun err, "\n" ]
+
+
+    --- Criterion ---
+
+    putStrLn "\n"
+
+
+    C.defaultMain
+       [ C.bench "com-error-1e-2" $ C.nf comFun 1e-2
+       , C.bench "com-error-1e-4" $ C.nf comFun 1e-4
+       , C.bench "com-error-1e-6" $ C.nf comFun 1e-6
+       , C.bench "com-error-1e-8" $ C.nf comFun 1e-8
+       , C.bench "com-error-1e-10" $ C.nf comFun 1e-10
+       , C.bench "com-error-1e-12" $ C.nf comFun 1e-12
+       , C.bench "com-error-1e-14" $ C.nf comFun 1e-14
+       , C.bench "expectation-step" $ C.nf (expectationStep zs) mx0
+       , C.bench "com-expectation-step" $ C.nf (expectationStep zs) cmx0
+       , C.bench "transition" $ C.nf toMean mx0
+       , C.bench "com-transition" $ C.nf toMean cmx0
+       , C.bench "mixture-step" $ C.nf (mixtureStep zs) mx0
+       , C.bench "com-mixture-step" $ C.nf (mixtureStep zs) cmx0 ]
diff --git a/goal-graphical.cabal b/goal-graphical.cabal
new file mode 100644
--- /dev/null
+++ b/goal-graphical.cabal
@@ -0,0 +1,74 @@
+cabal-version: 3.0
+version: 0.20
+name: goal-graphical
+synopsis: Optimization of latent variable and dynamical models with Goal
+description: goal-graphical provides tools for with dynamical and graphical models. Various graphical models are defined here, e.g. mixture models and restricted Boltzmann machines, dynamical models such as HMMs and Kalman filters, and in both cases algorithms for fitting them e.g. expectation maximization and contrastive divergence minimization.
+license: BSD-3-Clause
+license-file: LICENSE
+extra-source-files: README.md
+author: Sacha Sokoloski
+maintainer: sacha.sokoloski@mailbox.org
+homepage: https://gitlab.com/sacha-sokoloski/goal
+category: Math
+build-type: Simple
+
+library
+    exposed-modules:
+        Goal.Graphical.Models
+        Goal.Graphical.Models.Dynamic
+        Goal.Graphical.Models.Harmonium
+        Goal.Graphical.Models.Harmonium.FactorAnalysis
+        Goal.Graphical.Learning
+        Goal.Graphical.Inference
+        Goal.Graphical
+    build-depends:
+        base >= 4.13 && < 4.15,
+        mwc-random,
+        mwc-probability,
+        hmatrix-special,
+        ghc-typelits-knownnat,
+        ghc-typelits-natnormalise,
+        goal-core,
+        parallel,
+        statistics,
+        vector,
+        hmatrix,
+        containers,
+        goal-geometry,
+        goal-probability
+    default-language: Haskell2010
+    default-extensions:
+        NoStarIsType,
+        GeneralizedNewtypeDeriving,
+        StandaloneDeriving,
+        ScopedTypeVariables,
+        ExplicitNamespaces,
+        TypeOperators,
+        KindSignatures,
+        DataKinds,
+        RankNTypes,
+        TypeFamilies,
+        FlexibleContexts,
+        MultiParamTypeClasses,
+        ConstraintKinds,
+        FlexibleInstances
+    ghc-options: -Wall -O2
+
+benchmark com-poisson
+    type: exitcode-stdio-1.0
+    main-is: com-poisson.hs
+    hs-source-dirs: benchmarks
+    build-depends:
+        base,
+        goal-core,
+        goal-geometry,
+        goal-probability,
+        goal-graphical,
+        bytestring,
+        cassava,
+        criterion
+    default-language: Haskell2010
+    ghc-options: -Wall -O2
+    ghc-options: -threaded
+
+