diff --git a/CHANGES.md b/CHANGES.md
--- a/CHANGES.md
+++ b/CHANGES.md
@@ -1,6 +1,10 @@
 Revision history for Haskell package learning-hmm
 ===
 
+## Version 0.2.0.0
+- Remove dependency on the 'logfloat' package
+- Performance improvements
+
 ## Version 0.1.1.0
 - Add function `init` for random initialization
 - Add function `simulate` for running a Markov process
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.1.1.1
+version:             0.2.0.0
 stability:           experimental
 
 synopsis:            Yet another library for hidden Markov models
@@ -27,14 +27,13 @@
   exposed-modules:   Learning.HMM
   other-modules:     Data.Random.Distribution.Categorical.Util
                    , Data.Random.Distribution.Simplex
-                   , Data.Random.Distribution.Uniform.Util
-                   , Data.Vector.Util
-                   , Data.Vector.Util.LinearAlgebra
+                   , Data.Vector.Generic.Util
+                   , Data.Vector.Generic.Util.LinearAlgebra
                    , Learning.HMM.Internal
   -- other-extensions:  
   build-depends:     base >=4.7 && <4.8
                    , containers
-                   , logfloat
+                   , deepseq
                    , random-fu
                    , random-source
                    , vector
diff --git a/src/Data/Random/Distribution/Categorical/Util.hs b/src/Data/Random/Distribution/Categorical/Util.hs
--- a/src/Data/Random/Distribution/Categorical/Util.hs
+++ b/src/Data/Random/Distribution/Categorical/Util.hs
@@ -2,11 +2,10 @@
 
 module Data.Random.Distribution.Categorical.Util () where
 
+import Data.Maybe (fromMaybe)
 import Data.Random.Distribution (PDF, pdf)
-import Data.Random.Distribution.Categorical (Categorical, toList, totalWeight)
+import Data.Random.Distribution.Categorical (Categorical, toList)
 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
+  pdf cat a = fromMaybe 0 (lookup a $ map swap $ toList cat)
diff --git a/src/Data/Random/Distribution/Simplex.hs b/src/Data/Random/Distribution/Simplex.hs
--- a/src/Data/Random/Distribution/Simplex.hs
+++ b/src/Data/Random/Distribution/Simplex.hs
@@ -11,7 +11,6 @@
     , fractionalStdSimplex
     ) where
 
-import Control.Applicative
 import Control.Monad
 import Data.List
 import Data.Random.RVar
diff --git a/src/Data/Random/Distribution/Uniform/Util.hs b/src/Data/Random/Distribution/Uniform/Util.hs
deleted file mode 100644
--- a/src/Data/Random/Distribution/Uniform/Util.hs
+++ /dev/null
@@ -1,17 +0,0 @@
-{-# LANGUAGE FlexibleContexts, FlexibleInstances, MultiParamTypeClasses #-}
-
-module Data.Random.Distribution.Uniform.Util () where
-
-import Control.Applicative ((<$>))
-import Data.Number.LogFloat (LogFloat, logFloat, fromLogFloat)
-import Data.Random.Distribution (Distribution)
-import Data.Random.Distribution.Uniform -- (StdUniform(..), Uniform(..), doubleUniform)
-import Data.Random (rvarT)
-import Data.Random.Source (getRandomDouble)
-
-instance Distribution Uniform LogFloat where
-  rvarT (Uniform a b) = do x <- doubleUniform (fromLogFloat a) (fromLogFloat b)
-                           return $ logFloat x
-
-instance Distribution StdUniform LogFloat where
-  rvarT _ = logFloat <$> getRandomDouble
diff --git a/src/Data/Vector/Generic/Util.hs b/src/Data/Vector/Generic/Util.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Vector/Generic/Util.hs
@@ -0,0 +1,19 @@
+-- | Miscellaneous utility functions for "Data.Vector"
+module Data.Vector.Generic.Util (
+    frequencies
+  ) where
+
+import Data.Map.Strict (Map)
+import qualified Data.Map.Strict as M (empty, insertWith)
+import Data.Vector.Generic (Vector, foldl')
+
+-- $setup
+-- >>> :module + Data.Vector
+
+-- | @frequencies xs@ returns a 'Map' from distinct items in @xs@ to
+-- the number of times they appear.
+--
+-- >>> frequencies $ fromList "bra bra bar"
+-- fromList [(' ',2),('a',3),('b',3),('r',3)]
+frequencies :: (Ord a, Vector v a, Num n) => v a -> Map a n
+frequencies = foldl' (\m k -> M.insertWith (+) k 1 m) M.empty
diff --git a/src/Data/Vector/Generic/Util/LinearAlgebra.hs b/src/Data/Vector/Generic/Util/LinearAlgebra.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Vector/Generic/Util/LinearAlgebra.hs
@@ -0,0 +1,141 @@
+-- | 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.Generic.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.Generic (
+    Vector, any, cons, convert, 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 v a) => v a -> v a -> v 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 v a) => v a -> v a -> v 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 v a) => v a -> v a -> v 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 v a, Vector w (v a)) => w (v a) -> w (v a) -> w (v a)
+{-# INLINE (#+#) #-}
+m #+# n = zipWith (>+>) m n
+
+-- | Scalar-vector product
+--
+-- >>> 2 .> fromList [1, 2 :: Integer]
+-- fromList [2,4]
+(.>) :: (Num a, Vector v a) => a -> v a -> v a
+{-# INLINE (.>) #-}
+s .> v = map (s *) v
+
+-- | Vector-scalar division
+--
+-- >>> fromList [1, 2 :: Double] >/ 2
+-- fromList [0.5,1.0]
+(>/) :: (Fractional a, Vector v a) => v a -> a -> v 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 v a, Vector w (v a)) => w (v a) -> a -> w (v a)
+{-# INLINE (#/) #-}
+m #/ s = map (>/ s) m
+
+-- | Dot product
+--
+-- >>> fromList [1, 2] <.> fromList [3, 4 :: Int]
+-- 11
+(<.>) :: (Num a, Vector v a) => v a -> v 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 v a, Vector w (v a), Vector w a) => w (v a) -> v a -> v a
+{-# INLINE (#.>) #-}
+m #.> v = convert $ 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 v a, Vector w (v a), Vector w a) => v a -> w (v a) -> v a
+{-# INLINE (<.#) #-}
+v <.# m | any null m = empty
+        | otherwise  = (v <.> convert (map head m)) `cons` (v <.# map tail m)
+
+-- | Matrix transpose
+--
+-- >>> transpose $ fromList [fromList "ab", fromList "cd"]
+-- fromList [fromList "ac",fromList "bd"]
+transpose :: (Vector v a, Vector w (v a), Vector w a) => w (v a) -> w (v a)
+{-# INLINE transpose #-}
+transpose m
+  | any null m = empty
+  | otherwise  = convert (map head m) `cons` transpose (map tail m)
diff --git a/src/Data/Vector/Util.hs b/src/Data/Vector/Util.hs
deleted file mode 100644
--- a/src/Data/Vector/Util.hs
+++ /dev/null
@@ -1,13 +0,0 @@
--- | 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
deleted file mode 100644
--- a/src/Data/Vector/Util/LinearAlgebra.hs
+++ /dev/null
@@ -1,139 +0,0 @@
--- | 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
--- a/src/Learning/HMM.hs
+++ b/src/Learning/HMM.hs
@@ -11,7 +11,9 @@
 
 import Prelude hiding (init)
 import Control.Applicative ((<$>))
-import Control.Arrow ((***), first)
+import Control.Arrow (first)
+import Data.List (elemIndex, genericLength)
+import Data.Maybe (fromJust)
 import Data.Random.Distribution (pdf, rvar)
 import Data.Random.Distribution.Categorical (Categorical)
 import qualified Data.Random.Distribution.Categorical as C (
@@ -20,15 +22,12 @@
 import Data.Random.Distribution.Categorical.Util ()
 import Data.Random.RVar (RVar)
 import Data.Random.Sample (sample)
-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 qualified Data.Vector as V ((!), elemIndex, fromList, map, toList)
+import qualified Data.Vector.Generic as G (convert)
+import qualified Data.Vector.Generic.Util.LinearAlgebra as G (transpose)
+import qualified Data.Vector.Unboxed as U (fromList, toList)
 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' or 'init'.
 data HMM s o = HMM { states  :: [s] -- ^ Hidden states
@@ -89,17 +88,22 @@
 -- | @init states outputs@ returns a random variable of the model with
 --   @states@ and @outputs@, wherein parameters are sampled from uniform
 --   distributions.
-init :: (Ord s, Ord o) => [s] -> [o] -> RVar (HMM s o)
-init ss os = do hmm' <- init' (V.fromList ss) (V.fromList os)
-                return $ fromHMM' hmm'
+init :: (Eq s, Eq o) => [s] -> [o] -> RVar (HMM s o)
+init ss os = fromHMM' ss os <$> init' (length ss) (length os)
 
 -- | @model \`withEmission\` xs@ returns a model in which the
 --   'emissionDist' is updated by using the observed outputs @xs@. The
 --   'emissionDist' is set to be normalized histograms each of which is
 --   calculated from a partial set of @xs@ for each state. The partition is
 --   based on the most likely state path obtained by the Viterbi algorithm.
-withEmission :: (Ord s, Ord o) => HMM s o -> [o] -> HMM s o
-withEmission model xs = fromHMM' $ withEmission' (toHMM' model) (V.fromList xs)
+withEmission :: (Eq s, Eq o) => HMM s o -> [o] -> HMM s o
+withEmission model xs = fromHMM' ss os $ withEmission' model' xs'
+  where
+    ss     = states model
+    os     = outputs model
+    os'    = V.fromList os
+    model' = toHMM' model
+    xs'    = U.fromList $ fromJust $ mapM (`V.elemIndex` os') xs
 
 -- | @viterbi model xs@ performs the Viterbi algorithm using the observed
 --   outputs @xs@, and returns the most likely state path and its log
@@ -108,30 +112,36 @@
 viterbi model xs =
   checkModelIn "viterbi" model `seq`
   checkDataIn "viterbi" model xs `seq`
-  (V.toList *** logFromLogFloat) $ viterbi' model' xs'
+  first (V.toList . V.map (ss V.!) . G.convert) $ viterbi' model' xs'
   where
+    ss     = V.fromList $ states model
+    os'    = V.fromList $ outputs model
     model' = toHMM' model
-    xs'    = V.fromList xs
+    xs'    = U.fromList $ fromJust $ mapM (`V.elemIndex` os') xs
 
--- | @baumWelch model xs@ performs the Baum-Welch algorithm using the
---   observed outputs @xs@, and iteratively returns a list of updated
---   models and their corresponding log likelihoods.
+-- | @baumWelch model xs@ iteratively performs the Baum-Welch algorithm
+--   using the observed outputs @xs@, and returns a list of updated models
+--   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'
+  map (first $ fromHMM' ss os) $ baumWelch' model' xs'
   where
+    ss     = states model
+    os     = outputs model
+    os'    = V.fromList os
     model' = toHMM' model
-    xs'    = V.fromList xs
+    xs'    = U.fromList $ fromJust $ mapM (`V.elemIndex` os') xs
 
 -- | @simulate model t@ generates a Markov process of length @t@ using the
 --   @model@, and returns its state path and observed outputs.
 simulate :: HMM s o -> Int -> RVar ([s], [o])
-simulate model step | step < 1  = return ([], [])
-                    | otherwise = do s0 <- sample $ rvar pi0
-                                     x0 <- sample $ rvar $ phi s0
-                                     unzip . ((s0, x0) :) <$> sim s0 (step - 1)
+simulate model step
+  | step < 1  = return ([], [])
+  | otherwise = do s0 <- sample $ rvar pi0
+                   x0 <- sample $ rvar $ phi s0
+                   unzip . ((s0, x0) :) <$> sim s0 (step - 1)
   where
     sim _ 0 = return []
     sim s t = do s' <- sample $ rvar $ w s
@@ -164,33 +174,31 @@
     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
-                    }
+fromHMM' :: (Eq s, Eq o) => [s] -> [o] -> HMM' -> HMM s o
+fromHMM' ss os hmm' = HMM { states           = ss
+                          , outputs          = os
+                          , initialStateDist = C.fromList pi0'
+                          , transitionDist   = \s -> case elemIndex s ss of
+                                                       Nothing -> C.fromList []
+                                                       Just i  -> C.fromList $ w' i
+                          , emissionDist     = \s -> case 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
+    phi = G.transpose $ emissionDistT' hmm'
+    pi0'   = zip (U.toList pi0) ss
+    w' i   = zip (U.toList $ w V.! i) ss
+    phi' i = zip (U.toList $ phi V.! i) os
 
 -- | Convert 'HMM' to 'HMM''. The 'initialStateDist'', 'transitionDist'',
 --   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'
+toHMM' :: (Eq s, Eq o) => HMM s o -> HMM'
+toHMM' hmm = HMM' { nStates'          = length ss
+                  , nOutputs'         = length os
+                  , initialStateDist' = U.fromList pi0'
                   , transitionDist'   = V.fromList w'
                   , emissionDistT'    = V.fromList phi'
                   }
@@ -200,9 +208,9 @@
     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]
+    pi0' = [pdf pi0 s | s <- ss]
+    w'   = [U.fromList [pdf (w s) s' | s' <- ss] | s <- ss]
+    phi' = [U.fromList [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
--- a/src/Learning/HMM/Internal.hs
+++ b/src/Learning/HMM/Internal.hs
@@ -1,7 +1,6 @@
 module Learning.HMM.Internal (
     HMM' (..)
-  , Likelihood
-  , Probability
+  , LogLikelihood
   , init'
   , withEmission'
   , viterbi'
@@ -12,188 +11,206 @@
   ) where
 
 import Control.Applicative ((<$>))
+import Control.DeepSeq (NFData, force, rnf)
 import Control.Monad (forM_, replicateM)
 import Control.Monad.ST (runST)
-import qualified Data.Map.Strict as M (empty, insertWith, findWithDefault)
-import Data.Number.LogFloat (LogFloat, logFloat)
+import qualified Data.Map.Strict as M (findWithDefault)
 import Data.Random.RVar (RVar)
 import Data.Random.Distribution.Simplex (stdSimplex)
-import Data.Random.Distribution.Uniform.Util ()
-import Data.Vector (Vector, (!))
 import qualified Data.Vector as V (
-    filter, foldl', foldl1', freeze, fromList, last, length, map, maximum
-  , maxIndex, replicate, sum, tail, zip, zipWith, zipWith3, zipWith4
+    Vector, (!), filter, foldl1', freeze, fromList, generate, map , tail
+  , zip, zipWith, zipWith3
   )
+import qualified Data.Vector.Generic as G (convert)
+import qualified Data.Vector.Generic.Util as G (frequencies)
+import Data.Vector.Generic.Util.LinearAlgebra (
+    (>+>), (>.>), (>/>), (#+#), (.>), (>/), (#.>), (<.#)
+  )
+import qualified Data.Vector.Generic.Util.LinearAlgebra as G (transpose)
 import qualified Data.Vector.Mutable as MV (new, read, write)
-import qualified Data.Vector.Util as V (unsafeElemIndex)
-import Data.Vector.Util.LinearAlgebra (
-    (>+>), (>.>), (>/>), (#+#), (.>), (>/), (#/), (<.>), (#.>), (<.#)
+import qualified Data.Vector.Unboxed as U (
+    Vector, (!), freeze, fromList, generate, length, map, maxIndex, maximum
+  , replicate, sum, tail, zip
   )
-import qualified Data.Vector.Util.LinearAlgebra as V (transpose)
+import qualified Data.Vector.Unboxed.Mutable as MU (new, read, write)
 
-type Likelihood  = LogFloat
-type Probability = LogFloat
+type LogLikelihood = Double
 
--- | More efficient data structure of the 'HMM' model. 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)
-                     }
+-- | More efficient data structure of the 'HMM' model. The 'states' and
+--   'outputs' in 'HMM' are represented by their indices. The
+--   'initialStateDist', 'transitionDist', and 'emissionDist' are
+--   represented by matrices. The 'emissionDistT'' is a transposed matrix
+--   in order to simplify the calculation.
+data HMM' = HMM' { nStates'          :: Int -- ^ Number of states
+                 , nOutputs'         :: Int -- ^ Number of outputs
+                 , initialStateDist' :: U.Vector Double
+                 , transitionDist'   :: V.Vector (U.Vector Double)
+                 , emissionDistT'    :: V.Vector (U.Vector Double)
+                 }
 
-init' :: Vector s -> Vector o -> RVar (HMM' s o)
-init' ss os = do
-  let n = V.length ss
-      m = V.length os
-  pi0 <- V.fromList <$> stdSimplex (n-1)
-  w   <- V.fromList <$> replicateM n (V.fromList <$> stdSimplex (n-1))
-  phi <- V.fromList <$> replicateM n (V.fromList <$> stdSimplex (m-1))
-  return HMM' { states'           = ss
-              , outputs'          = os
+instance NFData HMM' where
+  rnf hmm' = rnf n `seq` rnf m `seq` rnf pi0 `seq` rnf w `seq` rnf phi'
+    where
+      n    = nStates' hmm'
+      m    = nOutputs' hmm'
+      pi0  = initialStateDist' hmm'
+      w    = transitionDist' hmm'
+      phi' = emissionDistT' hmm'
+
+init' :: Int -> Int -> RVar HMM'
+init' n m = do
+  pi0 <- U.fromList <$> stdSimplex (n-1)
+  w   <- V.fromList <$> replicateM n (U.fromList <$> stdSimplex (n-1))
+  phi <- V.fromList <$> replicateM n (U.fromList <$> stdSimplex (m-1))
+  return HMM' { nStates'          = n
+              , nOutputs'         = m
               , initialStateDist' = pi0
               , transitionDist'   = w
-              , emissionDistT'    = V.transpose phi
+              , emissionDistT'    = G.transpose phi
               }
 
-withEmission' :: (Ord s, Ord o) => HMM' s o -> Vector o -> HMM' s o
+withEmission' :: HMM' -> U.Vector Int -> HMM'
 withEmission' model xs = model { emissionDistT' = phi' }
   where
-    ss = states' model
-    os = outputs' model
-    (path, _) = viterbi' model xs
-    mp    = V.foldl' (\m k -> M.insertWith (+) k 1 m) M.empty $ V.zip path xs
-    hists = V.map (\s -> V.map (\o -> M.findWithDefault 0 (s, o) mp) os) ss
-    phi'  = V.transpose $ V.map (\h -> h >/ V.sum h) hists
+    ss   = V.generate (nStates' model) id
+    os   = U.generate (nOutputs' model) id
+    phi' = let (path, _) = viterbi' model xs
+               freqs     = G.frequencies $ U.zip path xs
+               hists     = V.map (\s -> U.map (\o ->
+                                 M.findWithDefault 0 (s, o) freqs) os) ss
+           in V.map (\f -> f >/ U.sum f) hists
 
-viterbi' :: Eq o => HMM' s o -> Vector o -> (Vector s, Likelihood)
-viterbi' model xs = (path, likelihood)
+viterbi' :: HMM' -> U.Vector Int -> (U.Vector Int, LogLikelihood)
+viterbi' model xs = (path, logL)
   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 <- MV.new n
-      ix `MV.write` (n-1) $ V.maxIndex $ deltas ! (n-1)
-      forM_ (reverse [0..(n-2)]) $ \i -> do
-        j <- ix `MV.read` (i+1)
-        ix `MV.write` i $ psis ! (i+1) ! j
-      V.freeze ix
-      where
-        ss = states' model
-    likelihood = V.maximum $ deltas ! (n-1)
+    n = U.length xs
 
-    deltas :: Vector (Vector Probability)
-    psis   :: Vector (Vector Int)
+    -- First, we calculate the value function and the state maximizing it
+    -- for each time.
+    deltas :: V.Vector (U.Vector Double)
+    psis   :: V.Vector (U.Vector Int)
     (deltas, psis) = runST $ do
       ds <- MV.new n
       ps <- MV.new n
-      ds `MV.write` 0 $ (phi' ! x 0) >.> pi0
-      ps `MV.write` 0 $ V.replicate k (0 :: Int)
-      forM_ [1..(n-1)] $ \i -> do
-        d <- ds `MV.read` (i-1)
-        let dws = V.map (d >.>) w'
-        ds `MV.write` i $ phi' ! x i >.> V.map V.maximum dws
-        ps `MV.write` i $ V.map V.maxIndex dws
+      MV.write ds 0 $ U.map log (phi' V.! (xs U.! 0)) >+> U.map log pi0
+      MV.write ps 0 $ U.replicate k 0
+      forM_ [1..(n-1)] $ \t -> do
+        d <- MV.read ds (t-1)
+        let dws = V.map (\wj -> d >+> U.map log wj) w'
+        MV.write ds t $ U.map log (phi' V.! (xs U.! t)) >+> G.convert (V.map U.maximum dws)
+        MV.write ps t $ G.convert (V.map U.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
+        k    = nStates' model
         pi0  = initialStateDist' model
-        w'   = V.transpose $ transitionDist' model
+        w'   = G.transpose $ transitionDist' model
         phi' = emissionDistT' model
 
-    -- Here we assumed that
-    n = V.length xs
+    -- The most likely path and corresponding log likelihood is as follows.
+    path = runST $ do
+      ix <- MU.new n
+      MU.write ix (n-1) $ U.maxIndex (deltas V.! (n-1))
+      forM_ (reverse [0..(n-2)]) $ \t -> do
+        i <- MU.read ix (t+1)
+        MU.write ix t $ psis V.! (t+1) U.! i
+      U.freeze ix
+    logL = U.maximum $ deltas V.! (n-1)
 
-baumWelch' :: (Eq s, Eq o) => HMM' s o -> Vector o -> [(HMM' s o, Likelihood)]
-baumWelch' model xs = zip ms $ tail ells
+baumWelch' :: HMM' -> U.Vector Int -> [(HMM', LogLikelihood)]
+baumWelch' model xs = zip models (tail logLs)
   where
-    (ms, ells) = unzip $ iterate ((`baumWelch1'` xs) . fst) (model, undefined)
+    n = U.length xs
+    step (m, _)     = baumWelch1' m n xs
+    (models, logLs) = unzip $ iterate step (model, undefined)
 
 -- | Perform one step of the Baum-Welch algorithm and return the updated
 --   model and the likelihood of the old model.
-baumWelch1' :: (Eq s, Eq o) => HMM' s o -> Vector o -> (HMM' s o, Likelihood)
-baumWelch1' model xs = (model', likelihood)
+baumWelch1' :: HMM' -> Int -> U.Vector Int -> (HMM', LogLikelihood)
+baumWelch1' model n xs = force (model', logL)
   where
-    model' = model { initialStateDist' = pi0
-                   , transitionDist'   = w
-                   , emissionDistT'    = phi'
-                   }
-    likelihood = V.last ells
-
-    -- First, we calculate the alpha and beta values using the
+    -- First, we calculate the alpha, beta, and scaling values using the
     -- forward-backward algorithm.
-    alphas = forward' model xs
-    betas  = backward' model xs
+    (alphas, cs) = forward' model n xs
+    betas        = backward' model n xs cs
 
-    -- 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 scaling values, we calculate the
+    -- posterior distribution, i.e., gamma and xi values.
+    (gammas, xis) = posterior' model n xs alphas betas cs
 
-    -- 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')
+    -- Using the gamma and xi values, we obtain the optimal initial state
+    -- probability vector, transition probability matrix, and emission
+    -- probability matrix.
+    pi0  = gammas V.! 0
+    w    = let ds = V.foldl1' (#+#) xis -- denominators
+               ns = V.map U.sum ds      -- numerators
+           in V.zipWith (>/) ds ns
+    phi' = let gs' o = V.map snd $ V.filter ((== o) . fst) $ V.zip (G.convert xs) gammas
+               ds    = V.foldl1' (>+>) . gs'  -- denominators
+               ns    = V.foldl1' (>+>) gammas -- numerators
+           in V.map (\o -> ds o >/> ns) os
       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
+        os = V.generate (nOutputs' model) id
 
-    -- Here we assumed that
-    os = outputs' model
+    -- We finally obtain the new model and the likelihood for the old model.
+    model' = model { initialStateDist' = pi0
+                   , transitionDist'   = w
+                   , emissionDistT'    = phi'
+                   }
+    logL = - (U.sum $ U.map log cs)
 
-forward' :: Eq o => HMM' s o -> Vector o -> Vector (Vector Probability)
-forward' model xs = runST $ do
-  v <- MV.new n
-  v `MV.write` 0 $ (phi' ! x 0) >.> pi0
-  forM_ [1..(n-1)] $ \i -> do
-    a <- v `MV.read` (i-1)
-    v `MV.write` i $ (phi' ! x i) >.> (a <.# w)
-  V.freeze v
+-- | Return alphas and scaling variables.
+forward' :: HMM' -> Int -> U.Vector Int -> (V.Vector (U.Vector Double), U.Vector Double)
+{-# INLINE forward' #-}
+forward' model n xs = runST $ do
+  as <- MV.new n
+  cs <- MU.new n
+  let a0 = (phi' V.! (xs U.! 0)) >.> pi0
+      c0 = 1 / U.sum a0
+  MV.write as 0 (c0 .> a0)
+  MU.write cs 0 c0
+  forM_ [1..(n-1)] $ \t -> do
+    a <- MV.read as (t-1)
+    let a' = (phi' V.! (xs U.! t)) >.> (a <.# w)
+        c' = 1 / U.sum a'
+    MV.write as t (c' .> a')
+    MU.write cs t c'
+  as' <- V.freeze as
+  cs' <- U.freeze cs
+  return (as', cs')
   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
 
-backward' :: Eq o => HMM' s o -> Vector o -> Vector (Vector Probability)
-backward' model xs = runST $ do
-  v <- MV.new n
-  v `MV.write` (n-1) $ V.replicate k $ logFloat (1 :: Double)
-  forM_ (reverse [0..(n-2)]) $ \i -> do
-    b <- v `MV.read` (i+1)
-    v `MV.write` i $ w #.> ((phi' ! x (i+1)) >.> b)
-  V.freeze v
+-- | Return betas using scaling variables.
+backward' :: HMM' -> Int -> U.Vector Int -> U.Vector Double -> V.Vector (U.Vector Double)
+{-# INLINE backward' #-}
+backward' model n xs cs = runST $ do
+  bs <- MV.new n
+  let bE = U.replicate k 1
+      cE = cs U.! (n-1)
+  MV.write bs (n-1) $ cE .> bE
+  forM_ (reverse [0..(n-2)]) $ \t -> do
+    b <- MV.read bs (t+1)
+    let b' = w #.> ((phi' V.! (xs U.! (t+1))) >.> b)
+        c' = cs U.! t
+    MV.write bs t $ c' .> b'
+  V.freeze bs
   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
+    k    = nStates' model
+    w    = transitionDist' model
+    phi' = emissionDistT' model
+
+-- | Return the posterior distribution.
+posterior' :: HMM' -> Int -> U.Vector Int -> V.Vector (U.Vector Double) -> V.Vector (U.Vector Double) -> U.Vector Double -> (V.Vector (U.Vector Double), V.Vector (V.Vector (U.Vector Double)))
+{-# INLINE posterior' #-}
+posterior' model _ xs alphas betas cs = (gammas, xis)
+  where
+    gammas = V.zipWith3 (\a b c -> a >.> b >/ c) alphas betas (G.convert cs)
+    xis    = V.zipWith3 (\a b x -> let w' = V.zipWith (.>) (G.convert a) w
+                                   in V.map ((phi' V.! x) >.> b >.>) w')
+               alphas (V.tail betas) (G.convert $ U.tail xs)
     w    = transitionDist' model
     phi' = emissionDistT' model
diff --git a/tests/doctests.hs b/tests/doctests.hs
--- a/tests/doctests.hs
+++ b/tests/doctests.hs
@@ -4,6 +4,6 @@
 
 main :: IO ()
 main = doctest [ "-isrc"
-               , "src/Data/Vector/Util/LinearAlgebra.hs"
+               , "src/Data/Vector/Generic/Util/LinearAlgebra.hs"
                , "src/Learning/HMM.hs"
                ]