diff --git a/hmm-hmatrix.cabal b/hmm-hmatrix.cabal
--- a/hmm-hmatrix.cabal
+++ b/hmm-hmatrix.cabal
@@ -1,5 +1,5 @@
 Name:                hmm-hmatrix
-Version:             0.0.1
+Version:             0.0.2
 Synopsis:            Hidden Markov Models using HMatrix primitives
 Description:
   Hidden Markov Models implemented using HMatrix data types and operations.
@@ -39,7 +39,7 @@
 Cabal-Version:       >=1.10
 
 Source-Repository this
-  Tag:         0.0.1
+  Tag:         0.0.2
   Type:        darcs
   Location:    http://hub.darcs.net/thielema/hmm-hmatrix
 
@@ -66,14 +66,15 @@
     hmatrix >=0.16 && <0.17,
     explicit-exception >=0.1.7 && <0.2,
     lazy-csv >=0.5 && <0.6,
-    random >=1.0 && <1.1,
-    transformers >= 0.2 && <0.5,
-    non-empty >=0.2.1 && <0.3,
-    semigroups >=0.8.4.1 && <0.17,
+    random >=1.0 && <1.2,
+    transformers >= 0.2 && <0.6,
+    non-empty >=0.2.1 && <0.4,
+    semigroups >=0.17 && <0.19,
     containers >=0.4.2 && <0.6,
     array >=0.4 && <0.6,
-    utility-ht >=0.0.10 && <0.1,
-    base >=4.5 && <4.8
+    utility-ht >=0.0.12 && <0.1,
+    deepseq >=1.3 && <1.5,
+    base >=4.5 && <5
   Hs-Source-Dirs:      src
   Default-Language:    Haskell2010
   GHC-Options:         -Wall
diff --git a/src/Math/HiddenMarkovModel.hs b/src/Math/HiddenMarkovModel.hs
--- a/src/Math/HiddenMarkovModel.hs
+++ b/src/Math/HiddenMarkovModel.hs
@@ -6,6 +6,7 @@
    Gaussian, GaussianTrained,
    uniform,
    generate,
+   probabilitySequence,
    Normalized.logLikelihood,
    Normalized.reveal,
 
@@ -45,6 +46,7 @@
 
 import qualified Data.NonEmpty as NonEmpty
 import qualified Data.Array as Array
+import Data.Traversable (Traversable, mapAccumL)
 import Data.Foldable (Foldable)
 import Data.Array (accumArray)
 
@@ -80,6 +82,19 @@
        }
 
 
+probabilitySequence ::
+   (Traversable f, Distr.EmissionProb distr,
+    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>
+   T distr prob -> f (State, emission) -> f prob
+probabilitySequence hmm =
+   snd
+   .
+   mapAccumL
+      (\index (State s, e) ->
+         (NC.atIndex (transition hmm) . flip (,) s,
+          index s * Distr.emissionStateProb (distribution hmm) e (State s)))
+      (NC.atIndex (initial hmm))
+
 generate ::
    (Rnd.RandomGen g, Ord prob, Rnd.Random prob,
     Distr.Generate distr, Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>
@@ -93,7 +108,7 @@
       return (x, takeColumn s $ transition hmm)
 
 takeColumn :: (Matrix.Element a) => Int -> Matrix a -> Vector a
-takeColumn n  =  Matrix.flatten . Matrix.extractRows [n] . Matrix.trans
+takeColumn n  =  Matrix.flatten . Matrix.extractColumns [n]
 
 
 
@@ -102,6 +117,7 @@
 -}
 trainSupervised ::
    (Distr.Estimate tdistr, Distr.Distribution tdistr ~ distr,
+    Distr.Trained distr ~ tdistr,
     Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>
    Int -> NonEmpty.T [] (State, emission) -> Trained tdistr prob
 trainSupervised n xs =
@@ -118,7 +134,7 @@
 
 finishTraining ::
    (Distr.Estimate tdistr, Distr.Distribution tdistr ~ distr,
-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>
+    Distr.Trained distr ~ tdistr, Distr.Probability distr ~ prob) =>
    Trained tdistr prob -> T distr prob
 finishTraining hmm =
    Cons {
@@ -131,7 +147,7 @@
 
 trainMany ::
    (Distr.Estimate tdistr, Distr.Distribution tdistr ~ distr,
-    Distr.Probability distr ~ prob,
+    Distr.Trained distr ~ tdistr, Distr.Probability distr ~ prob,
     Foldable f) =>
    (trainingData -> Trained tdistr prob) ->
    NonEmpty.T f trainingData -> T distr prob
diff --git a/src/Math/HiddenMarkovModel/CSV.hs b/src/Math/HiddenMarkovModel/CSV.hs
--- a/src/Math/HiddenMarkovModel/CSV.hs
+++ b/src/Math/HiddenMarkovModel/CSV.hs
@@ -16,6 +16,7 @@
 import Control.Monad.Exception.Synchronous (Exceptional)
 import Control.Monad (liftM2, replicateM, unless)
 
+import qualified Data.List.Reverse.StrictElement as Rev
 import qualified Data.List.HT as ListHT
 
 
@@ -113,7 +114,7 @@
    CSV.CSVRow -> CSVParser (Vector a)
 parseVectorFields =
    MT.lift . fmap Vector.fromList . mapM parseNumberCell .
-   ListHT.dropWhileRev (null . CSV.csvFieldContent)
+   Rev.dropWhile (null . CSV.csvFieldContent)
 
 parseNonEmptyVectorCells ::
    (Read a, Algo.Field a) =>
@@ -148,4 +149,4 @@
    CSV.CSVRow -> CSVParser [String]
 parseStringList =
    MT.lift . mapM cellContent .
-   ListHT.dropWhileRev (null . CSV.csvFieldContent)
+   Rev.dropWhile (null . CSV.csvFieldContent)
diff --git a/src/Math/HiddenMarkovModel/Distribution.hs b/src/Math/HiddenMarkovModel/Distribution.hs
--- a/src/Math/HiddenMarkovModel/Distribution.hs
+++ b/src/Math/HiddenMarkovModel/Distribution.hs
@@ -22,6 +22,7 @@
 import Numeric.Container ((<>))
 import Data.Packed.Matrix (Matrix)
 import Data.Packed.Vector (Vector)
+import Foreign.Storable (Storable)
 
 import qualified System.Random as Rnd
 
@@ -32,6 +33,7 @@
 import qualified Control.Monad.Exception.Synchronous as ME
 import qualified Control.Monad.Trans.Class as MT
 import qualified Control.Monad.Trans.State as MS
+import Control.DeepSeq (NFData, rnf)
 import Control.Monad (liftM2)
 
 import qualified Data.NonEmpty as NonEmpty
@@ -54,7 +56,10 @@
    toEnum = State
    fromEnum (State n) = n
 
+instance NFData State where
+   rnf (State n) = rnf n
 
+
 type family Probability distr
 type family Emission distr
 type family Trained distr
@@ -75,7 +80,13 @@
 class
    (NC.Container Vector (Probability distr), NC.Product (Probability distr)) =>
       EmissionProb distr where
+   {-
+   This function could be implemented generically in terms of emissionStateProb
+   but that would require an Info constraint.
+   -}
    emissionProb :: distr -> Emission distr -> Vector (Probability distr)
+   emissionStateProb :: distr -> Emission distr -> State -> Probability distr
+   emissionStateProb distr e (State s) = NC.atIndex (emissionProb distr e) s
 
 class
    (EmissionProb (Distribution tdistr),
@@ -102,7 +113,16 @@
 
 type instance Trained (Discrete prob symbol) = DiscreteTrained prob symbol
 
+
+instance (NFData prob, NFData symbol) => NFData (Discrete prob symbol) where
+   rnf (Discrete m) = rnf m
+
 instance
+   (NFData prob, NFData symbol) =>
+      NFData (DiscreteTrained prob symbol) where
+   rnf (DiscreteTrained m) = rnf m
+
+instance
    (NC.Container Vector prob, NC.Product prob, Ord symbol) =>
       Info (Discrete prob symbol) where
    numberOfStates (Discrete m) = Vector.dim $ snd $ Map.findMin m
@@ -160,10 +180,17 @@
 
 type instance Trained (Gaussian a) = GaussianTrained a
 
+
+instance (NFData a, Storable a) => NFData (Gaussian a) where
+   rnf (Gaussian params) = rnf params
+
+instance (NFData a, Storable a) => NFData (GaussianTrained a) where
+   rnf (GaussianTrained params) = rnf params
+
 instance (Algo.Field a) => Info (Gaussian a) where
    numberOfStates (Gaussian params) = Array.rangeSize $ Array.bounds params
 
-instance (Algo.Field a, Ord a, Rnd.Random a) => Generate (Gaussian a) where
+instance (Algo.Field a) => Generate (Gaussian a) where
    generate (Gaussian allParams) state = do
       let (center, covarianceChol, _c) = allParams ! state
       seed <- MS.state Rnd.random
@@ -174,13 +201,20 @@
             <> covarianceChol
 
 instance (HMatrix.Numeric a, Algo.Field a) => EmissionProb (Gaussian a) where
-   emissionProb (Gaussian allParams) =
-      let cholSolve m x =
-             Matrix.flatten $ Algo.cholSolve m $ Matrix.asColumn x
-          prob x (center, covarianceChol, c) =
-             let x0 = NC.sub x center
-             in  c * exp ((-1/2) * NC.dot x0 (cholSolve covarianceChol x0))
-      in  \x -> Vector.fromList $ map (prob x) $ Array.elems allParams
+   emissionProb (Gaussian allParams) x =
+      Vector.fromList $ map (emissionProbGen x) $ Array.elems allParams
+   emissionStateProb (Gaussian allParams) x s =
+      emissionProbGen x $ allParams ! s
+
+emissionProbGen ::
+   (HMatrix.Numeric a, Algo.Field a) =>
+   Vector a -> (Vector a, Matrix a, a) -> a
+emissionProbGen =
+   let cholSolve m x = Matrix.flatten $ Algo.cholSolve m $ Matrix.asColumn x
+   in  \x (center, covarianceChol, c) ->
+         let x0 = NC.sub x center
+         in  c * exp ((-1/2) * NC.dot x0 (cholSolve covarianceChol x0))
+
 
 instance (HMatrix.Numeric a, Algo.Field a) => Estimate (GaussianTrained a) where
    type Distribution (GaussianTrained a) = Gaussian a
diff --git a/src/Math/HiddenMarkovModel/Named.hs b/src/Math/HiddenMarkovModel/Named.hs
--- a/src/Math/HiddenMarkovModel/Named.hs
+++ b/src/Math/HiddenMarkovModel/Named.hs
@@ -21,6 +21,8 @@
 
 import qualified Control.Monad.Exception.Synchronous as ME
 import qualified Control.Monad.Trans.State as MS
+import Control.DeepSeq (NFData, rnf)
+import Foreign.Storable (Storable)
 
 import qualified Data.Map as Map
 import qualified Data.List as List
@@ -44,6 +46,12 @@
 
 type Discrete prob symbol = T (Distr.Discrete prob symbol) prob
 type Gaussian a = T (Distr.Gaussian a) a
+
+
+instance
+   (NFData distr, NFData prob, Storable prob) =>
+      NFData (T distr prob) where
+   rnf hmm = rnf (model hmm, nameFromStateMap hmm, stateFromNameMap hmm)
 
 
 fromModelAndNames :: HMM.T distr prob -> [String] -> T distr prob
diff --git a/src/Math/HiddenMarkovModel/Normalized.hs b/src/Math/HiddenMarkovModel/Normalized.hs
--- a/src/Math/HiddenMarkovModel/Normalized.hs
+++ b/src/Math/HiddenMarkovModel/Normalized.hs
@@ -59,7 +59,7 @@
 beta ::
    (Distr.EmissionProb distr,
     Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,
-    Traversable f, NonEmptyC.Zip f, NonEmptyC.Reverse f) =>
+    Traversable f, NonEmptyC.Reverse f) =>
    T distr prob ->
    f (prob, emission) -> NonEmpty.T f (Vector prob)
 beta hmm =
@@ -152,6 +152,7 @@
 -}
 trainUnsupervised ::
    (Distr.Estimate tdistr, Distr.Distribution tdistr ~ distr,
+    Distr.Trained distr ~ tdistr,
     Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>
    T distr prob -> NonEmpty.T [] emission -> Trained tdistr prob
 trainUnsupervised hmm xs =
diff --git a/src/Math/HiddenMarkovModel/Pattern.hs b/src/Math/HiddenMarkovModel/Pattern.hs
--- a/src/Math/HiddenMarkovModel/Pattern.hs
+++ b/src/Math/HiddenMarkovModel/Pattern.hs
@@ -37,7 +37,7 @@
 import Data.Packed.Vector (Vector)
 
 import qualified Data.Map as Map
-import Data.Semigroup (Semigroup, (<>), times1p)
+import Data.Semigroup (Semigroup, (<>), stimes)
 
 import Prelude hiding (replicate)
 
@@ -52,7 +52,7 @@
 
 instance (Algo.Field prob) => Semigroup (T prob) where
    (<>) = append
-   times1p k = replicate $ fromIntegral (k-1)
+   stimes k = replicate $ fromIntegral k
 
 
 infixl 5 `append`
diff --git a/src/Math/HiddenMarkovModel/Private.hs b/src/Math/HiddenMarkovModel/Private.hs
--- a/src/Math/HiddenMarkovModel/Private.hs
+++ b/src/Math/HiddenMarkovModel/Private.hs
@@ -17,6 +17,9 @@
 import Data.Packed.Matrix (Matrix)
 import Data.Packed.Vector (Vector)
 
+import Control.DeepSeq (NFData, rnf)
+import Foreign.Storable (Storable)
+
 import qualified Data.NonEmpty.Class as NonEmptyC
 import qualified Data.NonEmpty as NonEmpty
 import qualified Data.Semigroup as Sg
@@ -54,7 +57,12 @@
    }
    deriving (Show, Read)
 
+instance
+   (NFData distr, NFData prob, Storable prob) =>
+      NFData (T distr prob) where
+   rnf hmm = rnf (initial hmm, transition hmm, distribution hmm)
 
+
 emission ::
    (Distr.EmissionProb distr,
     Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>
@@ -140,8 +148,7 @@
       (NonEmpty.tail betas)
 
 zetaFromXi ::
-   (Distr.Probability distr ~ prob, Num prob, NC.Product prob) =>
-   T distr prob -> [Matrix prob] -> [Vector prob]
+   (NC.Product prob) => T distr prob -> [Matrix prob] -> [Vector prob]
 zetaFromXi hmm xis =
    map (NC.constant 1 (Matrix.rows $ transition hmm) <>) xis
 
@@ -212,18 +219,26 @@
    }
    deriving (Show, Read)
 
+instance
+   (NFData distr, NFData prob, Storable prob) =>
+      NFData (Trained distr prob) where
+   rnf hmm =
+      rnf (trainedInitial hmm, trainedTransition hmm, trainedDistribution hmm)
 
+
 sumTransitions ::
-   (NC.Container Vector e, Num e) =>
+   (NC.Container Vector e) =>
    T distr e -> [Matrix e] -> Matrix e
 sumTransitions hmm =
    List.foldl' NC.add (NC.konst 0 $ LinAlg.size $ transition hmm)
+--    zero = uncurry LinAlg.zeros $ LinAlg.size $ transition hmm
 
 {- |
 Baum-Welch algorithm
 -}
 trainUnsupervised ::
    (Distr.Estimate tdistr, Distr.Distribution tdistr ~ distr,
+    Distr.Trained distr ~ tdistr,
     Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>
    T distr prob -> NonEmpty.T [] emission -> Trained tdistr prob
 trainUnsupervised hmm xs =
@@ -243,7 +258,7 @@
 
 mergeTrained ::
    (Distr.Estimate tdistr, Distr.Distribution tdistr ~ distr,
-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>
+    Distr.Trained distr ~ tdistr, Distr.Probability distr ~ prob) =>
    Trained tdistr prob -> Trained tdistr prob -> Trained tdistr prob
 mergeTrained hmm0 hmm1 =
    Trained {
diff --git a/src/Math/HiddenMarkovModel/Test.hs b/src/Math/HiddenMarkovModel/Test.hs
--- a/src/Math/HiddenMarkovModel/Test.hs
+++ b/src/Math/HiddenMarkovModel/Test.hs
@@ -15,7 +15,9 @@
 
 import qualified Data.NonEmpty.Class as NonEmptyC
 import qualified Data.NonEmpty as NonEmpty
+import qualified Data.Traversable as Trav
 import qualified Data.Foldable as Fold
+import qualified Data.List as List
 import qualified Data.Map as Map
 import Data.NonEmpty ((!:))
 
@@ -43,13 +45,29 @@
 sequ :: NonEmpty.T [] Char
 sequ = 'a' !: take 20 (HMM.generate hmm (Rnd.mkStdGen 42))
 
+possibleStates :: Char -> [Distr.State]
+possibleStates c =
+   map Distr.State $ List.findIndices id $
+   map
+      (\p ->
+         case p of
+            0 -> False
+            1 -> True
+            _ -> error "invalid emission probability (must be 0 or 1)") $
+   Vector.toList $
+   Map.findWithDefault (error "invalid character") c $
+   case HMM.distribution hmm of Distr.Discrete m -> m
+
 {- |
 Should all be equal.
 -}
-sequLikelihood :: ((Double, Double), Double, NonEmpty.T [] Double)
+sequLikelihood :: ((Double, Double), Double, Double, NonEmpty.T [] Double)
 sequLikelihood =
    ((Priv.forward hmm sequ, Priv.backward hmm sequ),
     exp $ Normalized.logLikelihood hmm sequ,
+    sum $
+       map (NonEmpty.product . HMM.probabilitySequence hmm) $
+       Trav.mapM (\c -> map (flip (,) c) $ possibleStates c) sequ,
     NonEmptyC.zipWith NC.dot
        (Priv.alpha hmm sequ)
        (Priv.beta hmm $ NonEmpty.tail sequ))
diff --git a/src/Math/HiddenMarkovModel/Utility.hs b/src/Math/HiddenMarkovModel/Utility.hs
--- a/src/Math/HiddenMarkovModel/Utility.hs
+++ b/src/Math/HiddenMarkovModel/Utility.hs
@@ -10,12 +10,11 @@
 
 
 normalizeProb ::
-   (NC.Container Vector a, Fractional a) => Vector a -> Vector a
+   (NC.Container Vector a) => Vector a -> Vector a
 normalizeProb = snd . normalizeFactor
 
 normalizeFactor ::
-   (NC.Container Vector a, Fractional a) =>
-   Vector a -> (a, Vector a)
+   (NC.Container Vector a) => Vector a -> (a, Vector a)
 normalizeFactor xs =
    let c = NC.sumElements xs
    in  (c, NC.scale (recip c) xs)