diff --git a/CHANGES.md b/CHANGES.md
new file mode 100644
--- /dev/null
+++ b/CHANGES.md
@@ -0,0 +1,9 @@
+Revision history for Haskell package learning-hmm
+===
+
+## Version 0.1.1.0
+- Add function `init` for random initialization
+- Add function `simulate` for running a Markov process
+
+## Version 0.1.0.0
+- Original version
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.0.0
+version:             0.1.1.0
 stability:           experimental
 
 synopsis:            Yet another library for hidden Markov models
@@ -17,6 +17,7 @@
 
 cabal-version:       >=1.10
 build-type:          Simple
+extra-source-files:  CHANGES.md
 
 source-repository head
   type:              git
@@ -25,13 +26,17 @@
 library
   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
                    , Learning.HMM.Internal
   -- other-extensions:  
   build-depends:     base >=4.7 && <4.8
+                   , containers
                    , logfloat
                    , random-fu
+                   , random-source
                    , vector
   hs-source-dirs:    src
   default-language:  Haskell2010
diff --git a/src/Data/Random/Distribution/Simplex.hs b/src/Data/Random/Distribution/Simplex.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Random/Distribution/Simplex.hs
@@ -0,0 +1,50 @@
+{-# LANGUAGE
+    MultiParamTypeClasses,
+    FlexibleContexts, FlexibleInstances,
+    UndecidableInstances, GADTs
+  #-}
+
+module Data.Random.Distribution.Simplex
+    ( StdSimplex(..)
+    , stdSimplex
+    , stdSimplexT
+    , fractionalStdSimplex
+    ) where
+
+import Control.Applicative
+import Control.Monad
+import Data.List
+import Data.Random.RVar
+import Data.Random.Distribution
+import Data.Random.Distribution.Uniform
+
+-- |Uniform distribution over a standard simplex.
+newtype StdSimplex as =
+    -- | @StdSimplex k@ constructs a standard simplex of dimension @k@
+    -- (standard /k/-simplex).
+    -- An element of the simplex represents a vector variable @as = (a_0,
+    -- a_1, ..., a_k)@. The elements of @as@ are more than or equal to @0@
+    -- and @sum as@ is always equal to @1@.
+    StdSimplex Int
+    deriving (Eq, Show)
+
+instance (Ord a, Fractional a, Distribution StdUniform a) => Distribution StdSimplex [a] where
+    rvar (StdSimplex k) = fractionalStdSimplex k
+
+-- |@stdSimplex k@ returns a random variable being uniformly distributed over
+-- a standard simplex of dimension @k@.
+stdSimplex :: Distribution StdSimplex [a] => Int -> RVar [a]
+stdSimplex k = rvar (StdSimplex k)
+
+stdSimplexT :: Distribution StdSimplex [a] => Int -> RVarT m [a]
+stdSimplexT k = rvarT (StdSimplex k)
+
+-- |An algorithm proposed by Rubinstein & Melamed (1998).
+-- See, /e.g./, S. Onn, I. Weissman.
+-- Generating uniform random vectors over a simplex with implications to
+-- the volume of a certain polytope and to multivariate extremes.
+-- /Ann Oper Res/ (2011) __189__:331-342.
+fractionalStdSimplex :: (Ord a, Fractional a, Distribution StdUniform a) => Int -> RVar [a]
+fractionalStdSimplex k = do us <- replicateM k stdUniform
+                            let us' = sort us ++ [1]
+                            return $ zipWith (-) us' (0 : us')
diff --git a/src/Data/Random/Distribution/Uniform/Util.hs b/src/Data/Random/Distribution/Uniform/Util.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Random/Distribution/Uniform/Util.hs
@@ -0,0 +1,17 @@
+{-# 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/Learning/HMM.hs b/src/Learning/HMM.hs
--- a/src/Learning/HMM.hs
+++ b/src/Learning/HMM.hs
@@ -2,17 +2,24 @@
     HMM (..)
   , LogLikelihood
   , new
+  , init
+  , withEmission
   , viterbi
   , baumWelch
+  , simulate
   ) where
 
+import Prelude hiding (init)
+import Control.Applicative ((<$>))
 import Control.Arrow ((***), first)
-import Data.Random.Distribution (pdf)
+import Data.Random.Distribution (pdf, rvar)
 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.Random.RVar (RVar)
+import Data.Random.Sample (sample)
 import Data.List (genericLength)
 import Data.Number.LogFloat (fromLogFloat, logFloat, logFromLogFloat)
 import Data.Vector ((!))
@@ -23,11 +30,11 @@
 type LogLikelihood = Double
 
 -- | Parameter set of the hidden Markov model. Direct use of the
---   constructor is not recommended. Instead, call 'new'.
+--   constructor is not recommended. Instead, call 'new' or 'init'.
 data HMM s o = HMM { states  :: [s] -- ^ Hidden states
-                   , outputs :: [o] -- ^ Observed outputs
+                   , outputs :: [o] -- ^ Outputs
                    , initialStateDist :: Categorical Double s
-                     -- ^ Categorical distribusion of initial states
+                     -- ^ Categorical distribution of initial states
                    , transitionDist :: s -> Categorical Double s
                      -- ^ Categorical distribution of next states
                      --   conditioned by the previous states
@@ -53,11 +60,14 @@
     w   = transitionDist hmm
     phi = emissionDist hmm
 
--- | Construct a 'HMM' from the given states and outputs. The
---   'initialStateDist' and 'emissionDist' are set to be uniform
+-- | @new states outputs@ returns a model from the @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 [1, 2 :: Int] ['C', 'D']
+--   HMM {states = [1,2], outputs = "CD", initialStateDist = fromList [(0.5,1),(0.5,2)], transitionDist = [(fromList [(0.75,1),(0.25,2)],1),(fromList [(0.25,1),(0.75,2)],2)], emissionDist = [(fromList [(0.5,'C'),(0.5,'D')],1),(fromList [(0.5,'C'),(0.5,'D')],2)]}
 new :: (Ord s, Ord o) => [s] -> [o] -> HMM s o
 new ss os = HMM { states           = ss
                 , outputs          = os
@@ -76,8 +86,24 @@
     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.
+-- | @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'
+
+-- | @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)
+
+-- | @viterbi model xs@ performs the Viterbi algorithm using the observed
+--   outputs @xs@, and returns 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`
@@ -87,9 +113,9 @@
     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 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 :: (Eq s, Eq o) => HMM s o -> [o] -> [(HMM s o, LogLikelihood)]
 baumWelch model xs =
   checkModelIn "baumWelch" model `seq`
@@ -99,7 +125,23 @@
     model' = toHMM' model
     xs'    = V.fromList xs
 
--- | Check if the 'HMM' is valid in the sense of whether the 'states' and
+-- | @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)
+  where
+    sim _ 0 = return []
+    sim s t = do s' <- sample $ rvar $ w s
+                 x' <- sample $ rvar $ phi s
+                 ((s', x') :) <$> sim s' (t - 1)
+    pi0 = initialStateDist model
+    w   = transitionDist model
+    phi = emissionDist model
+
+-- | Check if the model is valid in the sense of whether the 'states' and
 --   'outputs' are not empty.
 checkModelIn :: String -> HMM s o -> ()
 checkModelIn fun hmm
@@ -111,8 +153,8 @@
     os = outputs hmm
     err = errorIn fun
 
--- | Check if all the elements of the data are contained in the 'outputs'
---   of the 'HMM'.
+-- | Check if all the elements of the observed outputs are contained in the
+--   'outputs' of the model.
 checkDataIn :: Eq o => String -> HMM s o -> [o] -> ()
 checkDataIn fun hmm xs
   | all (`elem` os) xs = ()
@@ -143,7 +185,7 @@
     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'',
+-- | 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
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
@@ -2,6 +2,8 @@
     HMM' (..)
   , Likelihood
   , Probability
+  , init'
+  , withEmission'
   , viterbi'
   , baumWelch'
   -- , baumWelch1'
@@ -9,15 +11,20 @@
   -- , backward'
   ) where
 
-import Control.Monad (forM_)
+import Control.Applicative ((<$>))
+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 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, foldl1', freeze, last, length, map, maximum, maxIndex
-  , replicate, sum , tail, zip , zipWith, zipWith3, zipWith4
+    filter, foldl', foldl1', freeze, fromList, 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.Mutable as MV (new, read, write)
 import qualified Data.Vector.Util as V (unsafeElemIndex)
 import Data.Vector.Util.LinearAlgebra (
     (>+>), (>.>), (>/>), (#+#), (.>), (>/), (#/), (<.>), (#.>), (<.#)
@@ -27,7 +34,7 @@
 type Likelihood  = LogFloat
 type Probability = LogFloat
 
--- | More efficient data structure of the HMM parameters. This should be
+-- | 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
@@ -37,19 +44,41 @@
                      , emissionDistT'    :: Vector (Vector Probability)
                      }
 
--- | Perform the Viterbi algorithm and return the most likely state path
---   and its likelihood.
+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
+              , initialStateDist' = pi0
+              , transitionDist'   = w
+              , emissionDistT'    = V.transpose phi
+              }
+
+withEmission' :: (Ord s, Ord o) => HMM' s o -> Vector o -> HMM' s o
+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
+
 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)
+      ix <- MV.new n
+      ix `MV.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
+        j <- ix `MV.read` (i+1)
+        ix `MV.write` i $ psis ! (i+1) ! j
       V.freeze ix
       where
         ss = states' model
@@ -58,15 +87,15 @@
     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)
+      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 `M.read` (i-1)
+        d <- ds `MV.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 `MV.write` i $ phi' ! x i >.> V.map V.maximum dws
+        ps `MV.write` i $ V.map V.maxIndex dws
       ds' <- V.freeze ds
       ps' <- V.freeze ps
       return (ds', ps')
@@ -82,15 +111,13 @@
     -- 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.
+--   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)
   where
@@ -137,14 +164,13 @@
     -- 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
+  v <- MV.new n
+  v `MV.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)
+    a <- v `MV.read` (i-1)
+    v `MV.write` i $ (phi' ! x i) >.> (a <.# w)
   V.freeze v
   where
     n   = V.length xs
@@ -155,14 +181,13 @@
     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)
+  v <- MV.new n
+  v `MV.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)
+    b <- v `MV.read` (i+1)
+    v `MV.write` i $ w #.> ((phi' ! x (i+1)) >.> b)
   V.freeze v
   where
     n   = V.length xs
diff --git a/tests/doctests.hs b/tests/doctests.hs
--- a/tests/doctests.hs
+++ b/tests/doctests.hs
@@ -1,6 +1,9 @@
+module Main (main) where
+
 import Test.DocTest
 
 main :: IO ()
 main = doctest [ "-isrc"
                , "src/Data/Vector/Util/LinearAlgebra.hs"
+               , "src/Learning/HMM.hs"
                ]