diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,22 @@
+The MIT License (MIT)
+
+Copyright (c) 2014 Mitsuhiro Nakamura
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
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/learning-hmm.cabal b/learning-hmm.cabal
new file mode 100644
--- /dev/null
+++ b/learning-hmm.cabal
@@ -0,0 +1,46 @@
+name:                learning-hmm
+version:             0.1.0.0
+stability:           experimental
+
+synopsis:            Yet another library for hidden Markov models
+description:         This library provides functions for the maximum likelihood
+                     estimation of discrete hidden Markov models. At present,
+                     only Baum-Welch and Viterbi algorithms are implemented.
+category:            Algorithms, Machine Learning, Statistics
+
+author:              Mitsuhiro Nakamura
+maintainer:          Mitsuhiro Nakamura <m.nacamura@gmail.com>
+copyright:           Copyright (c) 2014 Mitsuhiro Nakamura
+license:             MIT
+license-file:        LICENSE
+homepage:            https://github.com/mnacamura/learning-hmm
+
+cabal-version:       >=1.10
+build-type:          Simple
+
+source-repository head
+  type:              git
+  location:          https://github.com/mnacamura/learning-hmm.git
+
+library
+  exposed-modules:   Learning.HMM
+  other-modules:     Data.Random.Distribution.Categorical.Util
+                   , Data.Vector.Util
+                   , Data.Vector.Util.LinearAlgebra
+                   , Learning.HMM.Internal
+  -- other-extensions:  
+  build-depends:     base >=4.7 && <4.8
+                   , logfloat
+                   , random-fu
+                   , vector
+  hs-source-dirs:    src
+  default-language:  Haskell2010
+  ghc-options:       -Wall
+
+test-suite doctests
+  main-is:           doctests.hs
+  type:              exitcode-stdio-1.0
+  build-depends:     base, doctest >= 0.9.11
+  hs-source-dirs:    tests
+  default-language:  Haskell2010
+  ghc-options:       -threaded -Wall
diff --git a/src/Data/Random/Distribution/Categorical/Util.hs b/src/Data/Random/Distribution/Categorical/Util.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Random/Distribution/Categorical/Util.hs
@@ -0,0 +1,12 @@
+{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses #-}
+
+module Data.Random.Distribution.Categorical.Util () where
+
+import Data.Random.Distribution (PDF, pdf)
+import Data.Random.Distribution.Categorical (Categorical, toList, totalWeight)
+import Data.Tuple (swap)
+
+instance Eq a => PDF (Categorical Double) a where
+  pdf cat a = case lookup a $ map swap $ toList cat of
+                Nothing -> 0
+                Just p  -> p / totalWeight cat
diff --git a/src/Data/Vector/Util.hs b/src/Data/Vector/Util.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Vector/Util.hs
@@ -0,0 +1,13 @@
+-- | Miscellaneous utility functions for "Data.Vector"
+module Data.Vector.Util (
+    unsafeElemIndex
+  ) where
+
+import Data.Maybe (fromJust)
+import Data.Vector (Vector, elemIndex)
+
+-- | Return the index of the first occurrence of the given element or throw
+--   an error if no such occurrence.
+{-# INLINE unsafeElemIndex #-}
+unsafeElemIndex :: Eq a => a -> Vector a -> Int
+unsafeElemIndex e = fromJust . elemIndex e
diff --git a/src/Data/Vector/Util/LinearAlgebra.hs b/src/Data/Vector/Util/LinearAlgebra.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Vector/Util/LinearAlgebra.hs
@@ -0,0 +1,139 @@
+-- | Operators commonly used in the basic linear algebra. Note that all the
+--   functions defined here do not check the dimension/length of
+--   vectors/matrices.
+module Data.Vector.Util.LinearAlgebra (
+  -- * Pairwise operators
+    (>+>)
+  -- , (>->)
+  , (>.>)
+  , (>/>)
+  , (#+#)
+  -- , (#-#)
+  -- , (#.#)
+  -- , (#/#)
+
+  -- * Scalar-vector/vector-scalar operators
+  -- , (+>)
+  -- , (->)
+  , (.>)
+  -- , (/>)
+  -- , (>+)
+  -- , (>-)
+  -- , (>.)
+  , (>/)
+
+  -- * Scalar-matrix/matrix-scalar operators
+  -- , (+#)
+  -- , (-#)
+  -- , (.#)
+  -- , (/#)
+  -- , (#+)
+  -- , (#-)
+  -- , (#.)
+  , (#/)
+
+  -- * Dot and matrix-vector/vector-matrix products
+  , (<.>)
+  , (#.>)
+  , (<.#)
+
+  -- * Unary operators
+  , transpose
+  ) where
+
+import Prelude hiding (any, head, map, null, sum, tail, zipWith)
+import Data.Vector (Vector, any, cons, empty, head, map, null, sum, tail, zipWith)
+
+-- $setup
+-- >>> :module + Data.Vector
+
+-- | Pairwise addition between two vectors
+--
+-- >>> fromList [1, 2] >+> fromList [3, 4 :: Int]
+-- fromList [4,6]
+(>+>) :: Num a => Vector a -> Vector a -> Vector a
+{-# INLINE (>+>) #-}
+u >+> v = zipWith (+) u v
+
+-- | Pairwise product between two vectors
+--
+-- >>> fromList [1, 2] >.> fromList [3, 4 :: Double]
+-- fromList [3.0,8.0]
+(>.>) :: Num a => Vector a -> Vector a -> Vector a
+{-# INLINE (>.>) #-}
+u >.> v = zipWith (*) u v
+
+-- | Pairwise division between two vectors
+--
+-- >>> fromList [1, 2] >/> fromList [3, 4 :: Double]
+-- fromList [0.3333333333333333,0.5]
+(>/>) :: Fractional a => Vector a -> Vector a -> Vector a
+{-# INLINE (>/>) #-}
+u >/> v = zipWith (/) u v
+
+-- | Pairwise addition between two matrices
+--
+-- >>> fromList [fromList [1, 2], fromList [3, 4]] #+# fromList [fromList [5, 6], fromList [7, 8 :: Int]]
+-- fromList [fromList [6,8],fromList [10,12]]
+(#+#) :: Num a => Vector (Vector a) -> Vector (Vector a) -> Vector (Vector a)
+{-# INLINE (#+#) #-}
+m #+# n = zipWith (>+>) m n
+
+-- | Scalar-vector product
+--
+-- >>> 2 .> fromList [1, 2 :: Integer]
+-- fromList [2,4]
+(.>) :: Num a => a -> Vector a -> Vector a
+{-# INLINE (.>) #-}
+s .> v = map (s *) v
+
+-- | Vector-scalar division
+--
+-- >>> fromList [1, 2 :: Double] >/ 2
+-- fromList [0.5,1.0]
+(>/) :: Fractional a => Vector a -> a -> Vector a
+{-# INLINE (>/) #-}
+v >/ s = map (/ s) v
+
+-- | Matrix-scalar division
+--
+-- >>> fromList [fromList [1, 2], fromList [3, 4 :: Double]] #/ 2
+-- fromList [fromList [0.5,1.0],fromList [1.5,2.0]]
+(#/) :: Fractional a => Vector (Vector a) -> a -> Vector (Vector a)
+{-# INLINE (#/) #-}
+m #/ s = map (>/ s) m
+
+-- | Dot product
+--
+-- >>> fromList [1, 2] <.> fromList [3, 4 :: Int]
+-- 11
+(<.>) :: Num a => Vector a -> Vector a -> a
+{-# INLINE (<.>) #-}
+u <.> v = sum $ u >.> v
+
+-- | Matrix-vector product
+--
+-- >>> fromList [fromList [1, 2], fromList [3, 4]] #.> fromList [1, 2 :: Double]
+-- fromList [5.0,11.0]
+(#.>) :: Num a => Vector (Vector a) -> Vector a -> Vector a
+{-# INLINE (#.>) #-}
+m #.> v = map (<.> v) m
+
+-- | Vector-matrix product
+--
+-- >>> fromList [1, 2 :: Double] <.# fromList [fromList [1, 2], fromList [3, 4]]
+-- fromList [7.0,10.0]
+(<.#) :: Num a => Vector a -> Vector (Vector a) -> Vector a
+{-# INLINE (<.#) #-}
+v <.# m | any null m = empty
+        | otherwise  = (v <.> map head m) `cons` (v <.# map tail m)
+
+-- | Matrix transpose
+--
+-- >>> transpose $ fromList [fromList "ab", fromList "cd"]
+-- fromList [fromList "ac",fromList "bd"]
+transpose :: Vector (Vector a) -> Vector (Vector a)
+{-# INLINE transpose #-}
+transpose m
+  | any null m = empty
+  | otherwise  = map head m `cons` transpose (map tail m)
diff --git a/src/Learning/HMM.hs b/src/Learning/HMM.hs
new file mode 100644
--- /dev/null
+++ b/src/Learning/HMM.hs
@@ -0,0 +1,166 @@
+module Learning.HMM (
+    HMM (..)
+  , LogLikelihood
+  , new
+  , viterbi
+  , baumWelch
+  ) where
+
+import Control.Arrow ((***), first)
+import Data.Random.Distribution (pdf)
+import Data.Random.Distribution.Categorical (Categorical)
+import qualified Data.Random.Distribution.Categorical as C (
+    fromList, fromWeightedList, normalizeCategoricalPs
+  )
+import Data.Random.Distribution.Categorical.Util ()
+import Data.List (genericLength)
+import Data.Number.LogFloat (fromLogFloat, logFloat, logFromLogFloat)
+import Data.Vector ((!))
+import qualified Data.Vector as V (elemIndex, fromList, map, toList, zip)
+import qualified Data.Vector.Util.LinearAlgebra as V (transpose)
+import Learning.HMM.Internal
+
+type LogLikelihood = Double
+
+-- | Parameter set of the hidden Markov model. Direct use of the
+--   constructor is not recommended. Instead, call 'new'.
+data HMM s o = HMM { states  :: [s] -- ^ Hidden states
+                   , outputs :: [o] -- ^ Observed outputs
+                   , initialStateDist :: Categorical Double s
+                     -- ^ Categorical distribusion of initial states
+                   , transitionDist :: s -> Categorical Double s
+                     -- ^ Categorical distribution of next states
+                     --   conditioned by the previous states
+                   , emissionDist :: s -> Categorical Double o
+                     -- ^ Categorical distribution of outputs conditioned
+                     --   by the hidden states
+                   }
+
+instance (Show s, Show o) => Show (HMM s o) where
+  show = showHMM
+
+showHMM :: (Show s, Show o) => HMM s o -> String
+showHMM hmm = "HMM {states = "           ++ show ss
+              ++ ", outputs = "          ++ show os
+              ++ ", initialStateDist = " ++ show pi0
+              ++ ", transitionDist = "   ++ show [(w s, s) | s <- ss]
+              ++ ", emissionDist = "     ++ show [(phi s, s) | s <- ss]
+              ++ "}"
+  where
+    ss  = states hmm
+    os  = outputs hmm
+    pi0 = initialStateDist hmm
+    w   = transitionDist hmm
+    phi = emissionDist hmm
+
+-- | Construct a 'HMM' from the given states and outputs. The
+--   'initialStateDist' and 'emissionDist' are set to be uniform
+--   distributions. The 'transitionDist' is specified as follows: with
+--   probability 1/2, move to the same state, otherwise, move to a random
+--   state (which might be the same state).
+new :: (Ord s, Ord o) => [s] -> [o] -> HMM s o
+new ss os = HMM { states           = ss
+                , outputs          = os
+                , initialStateDist = pi0
+                , transitionDist   = w
+                , emissionDist     = phi
+                }
+  where
+    pi0 = C.fromWeightedList [(1, s) | s <- ss]
+    w s | s `elem` ss = C.fromList [(p s', s') | s' <- ss]
+        | otherwise   = C.fromList []
+      where
+        k = genericLength ss
+        p s' | s' == s   = 1/2 * (1 + 1/k)
+             | otherwise = 1/2 / k
+    phi s | s `elem` ss = C.fromWeightedList [(1, o) | o <- os]
+          | otherwise   = C.fromList []
+
+-- | Perform the Viterbi algorithm and return the most likely state path
+--   and its log likelihood.
+viterbi :: (Eq s, Eq o) => HMM s o -> [o] -> ([s], LogLikelihood)
+viterbi model xs =
+  checkModelIn "viterbi" model `seq`
+  checkDataIn "viterbi" model xs `seq`
+  (V.toList *** logFromLogFloat) $ viterbi' model' xs'
+  where
+    model' = toHMM' model
+    xs'    = V.fromList xs
+
+-- | Perform the Baum-Welch algorithm steps iteratively and return
+--   a list of updated 'HMM' parameters and their corresponding log
+--   likelihoods.
+baumWelch :: (Eq s, Eq o) => HMM s o -> [o] -> [(HMM s o, LogLikelihood)]
+baumWelch model xs =
+  checkModelIn "baumWelch" model `seq`
+  checkDataIn "baumWelch" model xs `seq`
+  map (fromHMM' *** logFromLogFloat) $ baumWelch' model' xs'
+  where
+    model' = toHMM' model
+    xs'    = V.fromList xs
+
+-- | Check if the 'HMM' is valid in the sense of whether the 'states' and
+--   'outputs' are not empty.
+checkModelIn :: String -> HMM s o -> ()
+checkModelIn fun hmm
+  | null ss   = err "empty states"
+  | null os   = err "empty outputs"
+  | otherwise = ()
+  where
+    ss = states hmm
+    os = outputs hmm
+    err = errorIn fun
+
+-- | Check if all the elements of the data are contained in the 'outputs'
+--   of the 'HMM'.
+checkDataIn :: Eq o => String -> HMM s o -> [o] -> ()
+checkDataIn fun hmm xs
+  | all (`elem` os) xs = ()
+  | otherwise          = err "illegal data"
+  where
+    os = outputs hmm
+    err = errorIn fun
+
+-- | Convert 'HMM'' to 'HMM'.
+fromHMM' :: (Eq s, Eq o) => HMM' s o -> HMM s o
+fromHMM' hmm' = HMM { states           = V.toList ss
+                    , outputs          = V.toList os
+                    , initialStateDist = C.fromList pi0'
+                    , transitionDist   = \s -> case V.elemIndex s ss of
+                                                 Nothing -> C.fromList []
+                                                 Just i  -> C.fromList $ w' i
+                    , emissionDist     = \s -> case V.elemIndex s ss of
+                                                 Nothing -> C.fromList []
+                                                 Just i  -> C.fromList $ phi' i
+                    }
+  where
+    ss  = states' hmm'
+    os  = outputs' hmm'
+    pi0 = initialStateDist' hmm'
+    w   = transitionDist' hmm'
+    phi = V.transpose $ emissionDistT' hmm'
+    pi0'   = V.toList $ V.map (first fromLogFloat) $ V.zip pi0 ss
+    w' i   = V.toList $ V.map (first fromLogFloat) $ V.zip (w ! i) ss
+    phi' i = V.toList $ V.map (first fromLogFloat) $ V.zip (phi ! i) os
+
+-- | Convert 'HMM' to 'HMM''. The 'initialStateDist'', 'transisionDist'',
+--   and 'emissionDistT'' are normalized.
+toHMM' :: (Eq s, Eq o) => HMM s o -> HMM' s o
+toHMM' hmm = HMM' { states'           = V.fromList ss
+                  , outputs'          = V.fromList os
+                  , initialStateDist' = V.fromList pi0'
+                  , transitionDist'   = V.fromList w'
+                  , emissionDistT'    = V.fromList phi'
+                  }
+  where
+    ss  = states hmm
+    os  = outputs hmm
+    pi0 = C.normalizeCategoricalPs $ initialStateDist hmm
+    w   = C.normalizeCategoricalPs . transitionDist hmm
+    phi = C.normalizeCategoricalPs . emissionDist hmm
+    pi0' = [logFloat $ pdf pi0 s | s <- ss]
+    w'   = [V.fromList [logFloat $ pdf (w s) s' | s' <- ss] | s <- ss]
+    phi' = [V.fromList [logFloat $ pdf (phi s) o | s <- ss] | o <- os]
+
+errorIn :: String -> String -> a
+errorIn fun msg = error $ "Learning.HMM." ++ fun ++ ": " ++ msg
diff --git a/src/Learning/HMM/Internal.hs b/src/Learning/HMM/Internal.hs
new file mode 100644
--- /dev/null
+++ b/src/Learning/HMM/Internal.hs
@@ -0,0 +1,174 @@
+module Learning.HMM.Internal (
+    HMM' (..)
+  , Likelihood
+  , Probability
+  , viterbi'
+  , baumWelch'
+  -- , baumWelch1'
+  -- , forward'
+  -- , backward'
+  ) where
+
+import Control.Monad (forM_)
+import Control.Monad.ST (runST)
+import Data.Number.LogFloat (LogFloat, logFloat)
+import Data.Vector (Vector, (!))
+import qualified Data.Vector as V (
+    filter, foldl1', freeze, last, length, map, maximum, maxIndex
+  , replicate, sum , tail, zip , zipWith, zipWith3, zipWith4
+  )
+import qualified Data.Vector.Mutable as M (new, read, write)
+import qualified Data.Vector.Util as V (unsafeElemIndex)
+import Data.Vector.Util.LinearAlgebra (
+    (>+>), (>.>), (>/>), (#+#), (.>), (>/), (#/), (<.>), (#.>), (<.#)
+  )
+import qualified Data.Vector.Util.LinearAlgebra as V (transpose)
+
+type Likelihood  = LogFloat
+type Probability = LogFloat
+
+-- | More efficient data structure of the HMM parameters. This should be
+--   only used internally. The 'emissionDistT'' is a transposed matrix in
+--   order to simplify the calculation.
+data HMM' s o = HMM' { states'           :: Vector s
+                     , outputs'          :: Vector o
+                     , initialStateDist' :: Vector Probability
+                     , transitionDist'   :: Vector (Vector Probability)
+                     , emissionDistT'    :: Vector (Vector Probability)
+                     }
+
+-- | Perform the Viterbi algorithm and return the most likely state path
+--   and its likelihood.
+viterbi' :: Eq o => HMM' s o -> Vector o -> (Vector s, Likelihood)
+viterbi' model xs = (path, likelihood)
+  where
+    -- The following procedure is based on
+    -- http://ibisforest.org/index.php?cmd=read&page=Viterbi%E3%82%A2%E3%83%AB%E3%82%B4%E3%83%AA%E3%82%BA%E3%83%A0&word=Viterbi
+    path = V.map (ss !) $ runST $ do
+      ix <- M.new n
+      ix `M.write` (n-1) $ V.maxIndex $ deltas ! (n-1)
+      forM_ (reverse [0..(n-2)]) $ \i -> do
+        j <- ix `M.read` (i+1)
+        ix `M.write` i $ psis ! (i+1) ! j
+      V.freeze ix
+      where
+        ss = states' model
+    likelihood = V.maximum $ deltas ! (n-1)
+
+    deltas :: Vector (Vector Probability)
+    psis   :: Vector (Vector Int)
+    (deltas, psis) = runST $ do
+      ds <- M.new n
+      ps <- M.new n
+      ds `M.write` 0 $ (phi' ! x 0) >.> pi0
+      ps `M.write` 0 $ V.replicate k (0 :: Int)
+      forM_ [1..(n-1)] $ \i -> do
+        d <- ds `M.read` (i-1)
+        let dws = V.map (d >.>) w'
+        ds `M.write` i $ phi' ! x i >.> V.map V.maximum dws
+        ps `M.write` i $ V.map V.maxIndex dws
+      ds' <- V.freeze ds
+      ps' <- V.freeze ps
+      return (ds', ps')
+      where
+        k   = V.length $ states' model
+        x i = let os  = outputs' model
+                  xs' = V.map (`V.unsafeElemIndex` os) xs
+              in xs' ! i
+        pi0  = initialStateDist' model
+        w'   = V.transpose $ transitionDist' model
+        phi' = emissionDistT' model
+
+    -- Here we assumed that
+    n = V.length xs
+
+-- | Perform the Baum-Welch algorithm steps iteratively and return
+--   a list of updated 'HMM'' parameters and their corresponding likelihoods.
+baumWelch' :: (Eq s, Eq o) => HMM' s o -> Vector o -> [(HMM' s o, Likelihood)]
+baumWelch' model xs = zip ms $ tail ells
+  where
+    (ms, ells) = unzip $ iterate ((`baumWelch1'` xs) . fst) (model, undefined)
+
+-- | Perform one step of the Baum-Welch algorithm and return the updated
+--   'HMM'' parameters and the likelihood of the old parameters.
+baumWelch1' :: (Eq s, Eq o) => HMM' s o -> Vector o -> (HMM' s o, Likelihood)
+baumWelch1' model xs = (model', likelihood)
+  where
+    model' = model { initialStateDist' = pi0
+                   , transitionDist'   = w
+                   , emissionDistT'    = phi'
+                   }
+    likelihood = V.last ells
+
+    -- First, we calculate the alpha and beta values using the
+    -- forward-backward algorithm.
+    alphas = forward' model xs
+    betas  = backward' model xs
+
+    -- Then, we obtain the likelihoods for each time. This should be
+    -- constant over time.
+    ells = V.zipWith (<.>) alphas betas
+
+    -- Based on the alpha, beta, and likelihood values, we calculate the
+    -- gamma and xi values.
+    gammas = V.zipWith3 (\a b l -> a >.> b >/ l) alphas betas ells
+    xis    = V.zipWith4 (\a b l x -> let w1 = V.zipWith (.>) a w0
+                                         w2 = V.map (phi0 ! x >.> b >.>) w1
+                                     in w2 #/ l)
+               alphas (V.tail betas) (V.tail ells) (V.tail xs')
+      where
+        xs'  = V.map (`V.unsafeElemIndex` os) xs
+        w0   = transitionDist' model
+        phi0 = emissionDistT' model
+
+    -- Using the gamma and xi values, we finally obtain the optimal initial
+    -- state probability vector, transition probability matrix, and
+    -- emission probability matrix.
+    pi0  = let gs = gammas ! 0
+           in gs >/ V.sum gs
+    w    = let ws = V.foldl1' (#+#) xis
+               zs = V.map V.sum ws
+           in V.zipWith (>/) ws zs
+    phi' = let gs' o = V.map snd $ V.filter ((== o) . fst) $ V.zip xs gammas
+               phis  = V.foldl1' (>+>) . gs'
+               zs    = V.foldl1' (>+>) gammas
+           in V.map (\o -> phis o >/> zs) os
+
+    -- Here we assumed that
+    os = outputs' model
+
+-- | Baum-Welch forward algorithm that generates α values
+forward' :: Eq o => HMM' s o -> Vector o -> Vector (Vector Probability)
+forward' model xs = runST $ do
+  v <- M.new n
+  v `M.write` 0 $ (phi' ! x 0) >.> pi0
+  forM_ [1..(n-1)] $ \i -> do
+    a <- v `M.read` (i-1)
+    v `M.write` i $ (phi' ! x i) >.> (a <.# w)
+  V.freeze v
+  where
+    n   = V.length xs
+    x i = let os  = outputs' model
+              xs' = V.map (`V.unsafeElemIndex` os) xs
+          in xs' ! i
+    pi0  = initialStateDist' model
+    w    = transitionDist' model
+    phi' = emissionDistT' model
+
+-- | Baum-Welch backward algorithm that generates β values
+backward' :: Eq o => HMM' s o -> Vector o -> Vector (Vector Probability)
+backward' model xs = runST $ do
+  v <- M.new n
+  v `M.write` (n-1) $ V.replicate k $ logFloat (1 :: Double)
+  forM_ (reverse [0..(n-2)]) $ \i -> do
+    b <- v `M.read` (i+1)
+    v `M.write` i $ w #.> ((phi' ! x (i+1)) >.> b)
+  V.freeze v
+  where
+    n   = V.length xs
+    k   = V.length $ states' model
+    x i = let os  = outputs' model
+              xs' = V.map (`V.unsafeElemIndex` os) xs
+          in xs' ! i
+    w    = transitionDist' model
+    phi' = emissionDistT' model
diff --git a/tests/doctests.hs b/tests/doctests.hs
new file mode 100644
--- /dev/null
+++ b/tests/doctests.hs
@@ -0,0 +1,6 @@
+import Test.DocTest
+
+main :: IO ()
+main = doctest [ "-isrc"
+               , "src/Data/Vector/Util/LinearAlgebra.hs"
+               ]