diff --git a/CHANGES.md b/CHANGES.md
--- a/CHANGES.md
+++ b/CHANGES.md
@@ -1,6 +1,9 @@
 Revision history for Haskell package learning-hmm
 ===
 
+## Version 0.3.2.1
+- Bug fix release
+
 ## Version 0.3.2.0
 - Add function `euclideanDistance` which measures the Euclidean distance between
   two models
diff --git a/learning-hmm.cabal b/learning-hmm.cabal
--- a/learning-hmm.cabal
+++ b/learning-hmm.cabal
@@ -1,5 +1,5 @@
 name:                learning-hmm
-version:             0.3.2.0
+version:             0.3.2.1
 stability:           experimental
 
 synopsis:            Yet another library for hidden Markov models
diff --git a/src/Learning/HMM.hs b/src/Learning/HMM.hs
--- a/src/Learning/HMM.hs
+++ b/src/Learning/HMM.hs
@@ -172,8 +172,8 @@
 --   'outputs' of the model.
 checkDataIn :: Eq o => String -> HMM s o -> [o] -> ()
 checkDataIn fun HMM {..} xs
-  | all (`elem` outputs) xs = ()
-  | otherwise               = errorIn fun "illegal data"
+  | any (`notElem` outputs) xs = errorIn fun "illegal data"
+  | otherwise                  = ()
 
 -- | Convert internal 'HMM' to 'HMM'.
 fromInternal :: (Eq s, Eq o) => [s] -> [o] -> I.HMM -> HMM s o
diff --git a/src/Learning/HMM/Internal.hs b/src/Learning/HMM/Internal.hs
--- a/src/Learning/HMM/Internal.hs
+++ b/src/Learning/HMM/Internal.hs
@@ -28,7 +28,7 @@
 import qualified Data.Vector.Mutable              as MV  ( unsafeNew, unsafeRead, unsafeWrite )
 import qualified Data.Vector.Unboxed              as U   ( Vector, fromList, length, map, sum, unsafeFreeze, unsafeIndex, unsafeTail, zip )
 import qualified Data.Vector.Unboxed.Mutable      as MU  ( unsafeNew, unsafeRead, unsafeWrite )
-import qualified Numeric.LinearAlgebra.Data       as H   ( (!), Matrix, Vector, diag, fromColumns, fromList, fromLists, fromRows, konst, maxElement, maxIndex, toColumns, tr )
+import qualified Numeric.LinearAlgebra.Data       as H   ( (!), Matrix, Vector, diag, fromColumns, fromList, fromLists, fromRows, ident, konst, maxElement, maxIndex, toColumns, tr )
 import qualified Numeric.LinearAlgebra.HMatrix    as H   ( (<>), (#>), sumElements )
 import           Prelude                          hiding ( init )
 
@@ -60,10 +60,12 @@
   phi <- H.fromLists <$> replicateM k (stdSimplex (l-1))
   return HMM { nStates          = k
              , nOutputs         = l
-             , initialStateDist = pi0
-             , transitionDist   = w
-             , emissionDistT    = H.tr phi
+             , initialStateDist = q_ H.#> pi0
+             , transitionDist   = w H.<> q_
+             , emissionDistT    = q_ H.<> H.tr phi
              }
+  where
+    q_ = q k -- Error matrix
 
 withEmission :: HMM -> U.Vector Int -> HMM
 withEmission (model @ HMM {..}) xs = model'
@@ -167,17 +169,34 @@
     -- posterior distribution, i.e., gamma and xi values.
     (gammas, xis) = posterior model n xs alphas betas cs
 
+    -- Error matrix
+    q_ = q nStates
+
     -- Using the gamma and xi values, we obtain the optimal initial state
     -- probability vector, transition probability matrix, and emission
     -- probability matrix.
-    pi0  = V.unsafeIndex gammas 0
+    pi0  = let g0  = V.unsafeIndex gammas 0
+               g0_ = g0 / H.konst (H.sumElements g0) nStates
+           in q_ H.#> g0_
     w    = let ds = V.foldl1' (+) xis         -- denominators
                ns = ds H.#> H.konst 1 nStates -- numerators
-           in H.diag (H.konst 1 nStates / ns) H.<> ds
+               w_ = H.diag (H.konst 1 nStates / ns) H.<> ds
+           in w_ H.<> q_
+           {- in H.fromRows $ zipWith3 (\n_ t t0 -> if n_ > eps then t else t0)
+            -                          (H.toList ns)
+            -                          (H.toRows $ w_ H.<> q_)
+            -                          (H.toRows transitionDist)
+            -}
     phi' = let gs' o = V.map snd $ V.filter ((== o) . fst) $ V.zip (G.convert xs) gammas
-               ds    = V.foldl' (+) 0 . gs'  -- denominators
+               ds    = V.foldl' (+) (H.konst 0 nStates) . gs'  -- denominators
                ns    = V.foldl1' (+) gammas -- numerators
-           in H.fromRows $ map (\o -> ds o / ns) [0..(nOutputs - 1)]
+               phi_  = H.fromRows $ map (\o -> ds o / ns) [0..(nOutputs - 1)]
+           in q_ H.<> phi_
+           {- in H.fromColumns $ zipWith3 (\n_ e e0 -> if n_ > eps then e else e0)
+            -                             (H.toList ns)
+            -                             (H.toColumns $ q_ H.<> phi_)
+            -                             (H.toColumns emissionDistT)
+            -}
 
     -- We finally obtain the new model and the likelihood for the old model.
     model' = model { initialStateDist = pi0
@@ -235,3 +254,26 @@
                alphas betas (G.convert cs)
     xis    = V.zipWith3 (\a b x -> H.diag a H.<> transitionDist H.<> H.diag (b * (emissionDistT H.! x)))
                alphas (V.unsafeTail betas) (G.convert $ U.unsafeTail xs)
+
+-- | Global error threshold.
+{-# INLINE eps #-}
+eps :: Double
+eps = 1e-4
+
+-- | Error matrix @q k@ is required to guarantee that the elements of initial
+--   states vector and emission/transition matrix are all larger than zero.
+--   @k@ is assumed to be the number of states. @q k@ is given by
+--       [    1 - eps, (1/k') eps,        ..., (1/k') eps ]
+--       [ (1/k') eps,    1 - eps,        ..., (1/k') eps ]
+--       [                    ...                         ]
+--       [ (1/k') eps,        ..., (1/k') eps,    1 - eps ],
+--   where the diagonal elements are @1 - eps@ and the remains are @(1/k')
+--   eps@. Here @eps@ is a small error value (given by @1e-4@) and
+--   @k' = k - 1@.
+q :: Int -> H.Matrix Double
+{-# INLINE q #-}
+q k = H.konst (1 - eps) (k, k) * e + H.konst (eps / k') (k, k) * (one - e)
+  where
+    e   = H.ident k
+    one = H.konst 1 (k, k)
+    k'  = fromIntegral (k - 1)
diff --git a/src/Learning/IOHMM.hs b/src/Learning/IOHMM.hs
--- a/src/Learning/IOHMM.hs
+++ b/src/Learning/IOHMM.hs
@@ -198,8 +198,10 @@
 --   in the 'inputs' ('outputs') of the model.
 checkDataIn :: (Eq i, Eq o) => String -> IOHMM i s o -> [i] -> [o] -> ()
 checkDataIn fun IOHMM {..} xs ys
-  | all (`elem` inputs) xs && all (`elem` outputs) ys = ()
-  | otherwise                                         = errorIn fun "illegal data"
+  | any (`notElem` inputs) xs  = errorIn fun "illegal input data"
+  | any (`notElem` outputs) ys = errorIn fun "illegal output data"
+  | any (`notElem` xs) inputs  = errorIn fun "insufficient input data"
+  | otherwise                  = ()
 
 -- | Convert internal 'IOHMM' to 'IOHMM'.
 fromInternal :: (Eq i, Eq s, Eq o) => [i] -> [s] -> [o] -> I.IOHMM -> IOHMM i s o
diff --git a/src/Learning/IOHMM/Internal.hs b/src/Learning/IOHMM/Internal.hs
--- a/src/Learning/IOHMM/Internal.hs
+++ b/src/Learning/IOHMM/Internal.hs
@@ -28,7 +28,7 @@
 import qualified Data.Vector.Mutable              as MV  ( unsafeNew, unsafeRead, unsafeWrite )
 import qualified Data.Vector.Unboxed              as U   ( Vector, fromList, length, map, sum, unsafeFreeze, unsafeIndex, unsafeTail, unzip, zip )
 import qualified Data.Vector.Unboxed.Mutable      as MU  ( unsafeNew, unsafeRead, unsafeWrite )
-import qualified Numeric.LinearAlgebra.Data       as H   ( (!), Matrix, Vector, diag, fromColumns, fromList, fromLists, fromRows, konst, maxElement, maxIndex, toColumns, tr )
+import qualified Numeric.LinearAlgebra.Data       as H   ( (!), Matrix, Vector, diag, fromColumns, fromList, fromLists, fromRows, ident, konst, maxElement, maxIndex, toColumns, tr )
 import qualified Numeric.LinearAlgebra.HMatrix    as H   ( (<>), (#>), sumElements )
 import           Prelude                          hiding ( init )
 
@@ -63,10 +63,12 @@
   return IOHMM { nInputs          = m
                , nStates          = k
                , nOutputs         = l
-               , initialStateDist = pi0
-               , transitionDist   = w
-               , emissionDistT    = H.tr phi
+               , initialStateDist = q_ H.#> pi0
+               , transitionDist   = V.map (H.<> q_) w
+               , emissionDistT    = q_ H.<> H.tr phi
                }
+  where
+    q_ = q k -- Error matrix
 
 withEmission :: IOHMM -> U.Vector (Int, Int) -> IOHMM
 withEmission (model @ IOHMM {..}) xys = model'
@@ -175,18 +177,37 @@
     -- posterior distribution, i.e., gamma and xi values.
     (gammas, xis) = posterior model n xys alphas betas cs
 
+    -- Error matrix
+    q_ = q nStates
+
     -- Using the gamma and xi values, we obtain the optimal initial state
     -- probability vector, transition probability matrix, and emission
     -- probability matrix.
-    pi0  = V.unsafeIndex gammas 0
+    -- Each simplex in pi0, w, and phi' remains old if their numerators are
+    -- zero.
+    pi0  = let g0  = V.unsafeIndex gammas 0
+               g0_ = g0 / H.konst (H.sumElements g0) nStates
+           in q_ H.#> g0_
     w    = let xis' i = V.map snd $ V.filter ((== i) . fst) $ V.zip (G.convert $ U.unsafeTail xs) xis
-               ds     = V.foldl1' (+) . xis'  -- denominators
+               ds     = V.foldl1' (+) . xis'        -- denominators
                ns i   = ds i H.#> H.konst 1 nStates -- numerators
-           in V.map (\i -> H.diag (H.konst 1 nStates / ns i) H.<> ds i) (V.generate nInputs id)
+               w_ i   = H.diag (H.konst 1 nStates / ns i) H.<> ds i
+           in flip V.map (V.generate nInputs id) $ \i -> w_ i H.<> q_
+                {- H.fromRows $ zipWith3 (\n_ t t0 -> if n_ > eps then t else t0)
+                 -                       (H.toList $ ns i)
+                 -                       (H.toRows $ w_ i H.<> q_)
+                 -                       (H.toRows $ V.unsafeIndex transitionDist i)
+                 -}
     phi' = let gs' o = V.map snd $ V.filter ((== o) . fst) $ V.zip (G.convert ys) gammas
-               ds    = V.foldl' (+) 0 . gs'  -- denominators
-               ns    = V.foldl1' (+) gammas -- numerators
-           in H.fromRows $ map (\o -> ds o / ns) [0..(nOutputs - 1)]
+               ds    = V.foldl' (+) (H.konst 0 nStates) . gs' -- denominators
+               ns    = V.foldl1' (+) gammas                   -- numerators
+               phi_  = H.fromRows $ map (\o -> ds o / ns) [0..(nOutputs - 1)]
+           in q_ H.<> phi_
+           {- in H.fromColumns $ zipWith3 (\n_ e e0 -> if n_ > eps then e else e0)
+            -                             (H.toList ns)
+            -                             (H.toColumns $ q_ H.<> phi_)
+            -                             (H.toColumns emissionDistT)
+            -}
 
     -- We finally obtain the new model and the likelihood for the old model.
     model' = model { initialStateDist = pi0
@@ -247,3 +268,26 @@
     xis    = V.zipWith3 (\a b (x, y) -> H.diag a H.<> w x H.<> H.diag (b * (emissionDistT H.! y)))
                alphas (V.unsafeTail betas) (G.convert $ U.unsafeTail xys)
     w = V.unsafeIndex transitionDist
+
+-- | Global error threshold.
+{-# INLINE eps #-}
+eps :: Double
+eps = 1e-4
+
+-- | Error matrix @q k@ is required to guarantee that the elements of initial
+--   states vector and emission/transition matrix are all larger than zero.
+--   @k@ is assumed to be the number of states. @q k@ is given by
+--       [    1 - eps, (1/k') eps,        ..., (1/k') eps ]
+--       [ (1/k') eps,    1 - eps,        ..., (1/k') eps ]
+--       [                    ...                         ]
+--       [ (1/k') eps,        ..., (1/k') eps,    1 - eps ],
+--   where the diagonal elements are @1 - eps@ and the remains are @(1/k')
+--   eps@. Here @eps@ is a small error value (given by @1e-4@) and
+--   @k' = k - 1@.
+q :: Int -> H.Matrix Double
+{-# INLINE q #-}
+q k = H.konst (1 - eps) (k, k) * e + H.konst (eps / k') (k, k) * (one - e)
+  where
+    e   = H.ident k
+    one = H.konst 1 (k, k)
+    k'  = fromIntegral (k - 1)