Hidden Markov model definitions

This notebook defines a number of functions related to hidden Markov models. There are two main visualization functions defined: (1) DrawHMM[ ] and (2) DrawObservationSequence; there are to main import functions defined: (1) ReadHMM[ ] and (2) ReadObservationSequences[ ] -- these functions assume the same format as 'hmm_train', 'hmm_eval' and 'hmm_viterbi' (see source code on web site); finally a number of HMM algorithmic functions are defined: (1) FloorHMM[ ], (2) RandomHMM[ ], (3) GenerateObservationSequence, (4) LogProb[ ], (5) NormProb[ ], (6) DecodeHMM[ ] and (7) TrainHMM (function name indicate the main role of each function).

General definitions

Turn off annoying warning messages

Off[General::spell];
Off[General::spell1];

Load packages

<< Graphics`Colors`;
<< Statistics`DataManipulation`;
<< Graphics`Graphics`;
<< Statistics`DescriptiveStatistics`;

Visualization definitions

Drawing subroutines

eq1  =  (x1  -  cx)^2  +  (y1  -  cy)^2  ==  
          r^2;
eq2  =  (x2  -  cx)^2  +  (y2  -  cy)^2  ==  
          r^2;
\(\(\(eq3\  = \ \((x3\  - \ cx)\)^2\  + \ \((y3\  - \ cy)\)^2\  == \ 
          r^2;\)\(
\)
    \)\)
\(\(\(sol\  = \ Flatten[Solve[{eq1, eq2, eq3}, \ {cx, cy, r}]];\)\(
\)
    \)\)
circleParameters[{{x1s_, y1s_},  {x2s_, y2s_},  {x3s_,  
              y3s_}}] := ({{cx,  cy},  Abs[r]}  /.  sol)  /.  {x1 -> 
              x1s,  x2 -> x2s, x3 -> x3s,  y1 -> y1s, y2 -> y2s, 
            y3 -> y3s} ;

eq4  = (x  -  cx1)^2  +  (y - cy1)^2  ==  
          r1^2;
\(\(\(eq5\  = \((x\  - \ cx2)\)^2\  + \ \((y - cy2)\)^2\  == \ 
          r2^2;\)\(
\)
    \)\)
\(\(\(sol2\  = \ Solve[{eq4, eq5}, \ {x, y}];\)\(
\)
    \)\)
circleIntersections[{{cx1s_, cy1s_}, r1s_}, {{cx2s_, cy2s_},  
            r2s_}] := ({x, y}  /.  sol2)  /.  {cx1 -> cx1s, 
            cx2 -> cx2s, cy1 -> cy1s, cy2 -> cy2s, r1 -> r1s,  
            r2 -> r2s};

DiskList[c_, r_, pr_,  sty_] := 
        Table[{sty[[i]],  Disk[c, r, Take[pr, {i, i + 1}]]},  {i,  
            1,  Length[pr] - 1}];

TransformVector[
          pr_] := (2  
              Pi ({0. }~Join~
                  Table[Plus  @@  Take[pr, {1, i}],  {i,  1,  
                      Length[pr]}])) //  N;

TextConvert[{{{cx_, cy_}, r_, {t1_, t2_}}, {i_, j_}},  
          offset_] := 
Module[{ t,  
            pos}, 
t  =  (t1  +  t2)/
              2; 
pos  =  ({cx, cy}  +  
                  offset r { Cos[t],  Sin[t]})  // N; 
          Text[Subscript["a", ToString[j] <> ToString[i]],  pos, 
            Background -> GrayLevel[1. ]]
];

ArrowHead[{α_,  β_,  
            h_},  {{x_, y_}, θ_}] := 
Module[{l1,  
            l2,  pts}, 
l1  =  
            h/(Cos[α]  +  
                  Sin[α] Tan[β]); 
l2  =  
            l1 Sin[α]/
                Cos[β]; 
pts  = {{0,  -l2 Sin[β]},  {l2 Cos[β],  0},  {0, 
                  l1 Cos[α]}, {-l2 Cos[β],  0}}  //  N; 
          pts  =  (# - {0, l1 Cos[α]}) &  /@  
              pts; ({{Cos[θ],  -Sin[θ]},  {Sin[θ],  Cos[θ]}} . #  +  {x, y}) &  /@ 
            pts
];

convert[{{{cx_, cy_}, r_}, θ_}] := {{cx  +  
              r Cos[θ],  
            cy  +  r Sin[θ]},  θ};

getOption[name_,  default_,  opts___] := 
        Module[{var}, 
          If[NumberQ[var  = (name  /.  {opts})  //  N],  var,  
            default]];

getOption2[name_,  default_,  opts___] := 
        Module[{var}, 
          If[Length[var  = (name  /.  {opts})]  >  0,  var,  
            default]];

getOption3[name_,  default_,  
          opts___] := 
Module[{var}, 
          If[(var  = (name  /.  {opts}))  ==  name,  default,  
            var, var]];

eliminateNulls[l_] := Delete[l,  Position[l,  Null]];

thirdPoint[{{x1_,  y1_}, {x2_, 
              y2_}},  α_] := 
Module[{dx,  
            dy,  θ,  
            b}, 
dx  =  (x2 - 
                x1); 
dy  =  (y2 - 
                y1); 
b  =  
            Sqrt[(x2 - x1)^2  +  (y2 - y1)^2]/
              2; 
θ  =  
            ArcTan[dx,  dy]  +  α; 
{{x1, 
                y1}, {x2, y2}, {x1  +  b Cos[θ],  
                y1  +  b Sin[θ]}}  //  
            N
];

arcAngles[{{cx_, cy_},  
            r_},  {{x1_,  y1_}, {x2_, 
              y2_}}] := 
Module[{dx,  dy,  t1,  t2, tt1, 
            tt2}, 
dx  =  (x1  -  
                cx); 
dy  =  (y1  -  
                cy); 
t1  =  
            ArcTan[dx,  
              dy]; 

dx  =  (x2  -  
                cx); 
dy  =  (y2  -  
                cy); 
t2  =  
            ArcTan[dx,  dy]; 

{t1, 
              t2}  =  Sort[{t1, t2}]; 
If[
            t2  -  t1  >  Pi,  tt1  = t2;  tt2  = t1 + 2  Pi, 
            tt1 = t1; tt2 = t2]; 
{{cx, cy}, 
            r, {tt1, tt2}}
];

\(\(selfAngles[{{cx_, cy_}, \ 
            r_}, \ {{x1_, \ y1_}, {x2_, 
              y2_}}] := 
Module[{dx, \ dy, \ t1, \ t2, tt1, 
            tt2}, 
dx\  = \ \((x1\  - \ 
                cx)\); 
dy\  = \ \((y1\  - \ 
                cy)\); 
t1\  = \ 
            ArcTan[dx, \ 
              dy]; 

dx\  = \ \((x2\  - \ 
                cx)\); 
dy\  = \ \((y2\  - \ 
                cy)\); 
t2\  = \ 
            ArcTan[dx, \ dy]; {t1, t2}\  = \ 
            Sort[{t1, t2}]; 

If[
            t2\  - \ t1\  < \ Pi, \ tt1\  = t2; \ tt2\  = t1 + 2  Pi, 
            tt1 = t1; tt2 = t2]; 
{{cx, cy}, 
            r, {tt1, tt2}}
];\)\)

minDistance[{{{xi1_, yi1_}, {xi2_, yi2_}}, {x3_, 
              y3_}}] := 
Module[{d1, 
            d2}, 
d1  =  
            Sqrt[(x3  -  xi1)^2  +  (y3  -  
                      yi1)^2]; 
d2  =  
            Sqrt[(x3  -  xi2)^2  +  (y3  -  
                      yi2)^2]; 
If[
            d1  <  d2,  {xi1,  yi1},  {xi2,  
              yi2}]
];

FixTextInternal[gr_,  
          l_] := 
Module[{gr2,  s1,  s2, 
            i}, 
gr2  =  
            gr; 
Do[{s1, s2}  =  
              l[[i]]; 
gr2  =  
              ReplacePart[gr2, s2, Position[gr2,  s1]],  {i,  1,  
              Length[l]}]; gr2
];

FixText[gr_,  l_] := 
Module[{gr2, t1, t2, t3}, 
          t1  =  {ToString[#[[1]]] <> ToString[#[[2]]], "t" <> 
                    ToString[#[[2]]] <> ToString[#[[1]]]} &  /@  
              l; 
gr2  =  FixTextInternal[gr, t1]; 
          t2  =  {ToString[#[[2]]] <> ToString[#[[1]]], 
                  ToString[#[[1]]] <> ToString[#[[2]]]} &  /@  l; 
          gr2  =  FixTextInternal[gr2, t2]; 
          t3  = {"t" <> ToString[#[[2]]] <> 
                    ToString[#[[1]]],  
                  ToString[#[[2]]] <> ToString[#[[1]]]} &  /@  
              l; 
Show[
            FixTextInternal[gr2,  t3]]
];

SubstituteValue[gr_, {i_, j_},  
          s_] := 
ReplacePart[gr,  s,  
          Position[gr,  
            Subscript["a", ToString[i] <> ToString[j]]]];

SubstituteValues[gr_,  
          a_] := 
Module[{gr2,  i, j, 
            s}, 
gr2  =  
            gr; 
Do[
s  =  
              ToString[
                InputForm[a[[i, j]]]]; 
gr2  =  
              SubstituteValue[gr2,  {i, j},  s],  {i,  1,  
              Length[a]},  {j,  1,  Length[a]}]; 
          Show[gr2]
];

Markov model drawing function (DrawMarkov)

Description

DrawMarkov[A, opts] takes the state transition matrix   A and zero or more options opt, and draws the corresponding Markov model. When the Markov model is labeled, some subscripts may be reversed. In that case,  FixText[Graphic, {{i1,j1},{i2,j2},...}] will switch the {ik,jk} subscripts. The function SubstituteValues[Graphic, A] will substitute actual values of the matrix A into the Markov model plot.

Options:

* MarkovStyle -> {Circular, Sequential}: arrange states  in a circle or in a line.
* AngularOffset -> 0 Degree:  offset for circularly arranged Markov models.
* StateSize -> 1/4: relative size of state circles.
* StateColor -> GrayLevel[.75]: color of each state.
* LineThickness -> 1/2 (point): default line thickness  for drawing.
* BendAngle -> 20 Degree:  curvature of state-transition arcs.
* SelfCircleSize -> 1/5: relative radius of self-connections.
* SelfCircleOffset -> 1/3: relative offset of self-connections from center of states.
* ArrowAngleTop -> 10 Degree: arrow tip spread angle.
* ArrowAngleBottom -> 10 Degree: arrow tail  angle.
* ArrowLength -> 0.075: relative length of arrows.
* TextOffset -> 1.1: transition arcs offset for text labeling.
* SelfTextOffset -> 0.1: self-transition arcs offset for text labeling.
* MarkovText -> True: whether to label Markov model with text.
* FontSize -> 12: what font size to use.
* FontFamily -> "Times": what font family to use.
* FontSlant -> "Italic": what font slant to use.

Definition

\(\(DrawMarkov[a_, \ 
          opts___] := 
Module[{at, c, \ n, \ α, \ 
            r, \ lt, \ g0, g1, g2, g3, g4, g5, g6, \ glist, lines, \ 
            linesTmp, \ h, \ arcPts, \ arcs, \ states, \ candidates, third, \ 
            intersect, arcsComplete, \ arrows, \ arrowAngle1, \ arrowAngle2, 
            arrowLength, \ selfCenters, \ r2, \ r3, \ intersect2, \ selfArcs, 
            selfComplete, \ thetas2, \ arrows2, \ stateColor, \ hmmStyle, \ 
            stateTxt, nonZero, \ nonZeroSelf, \ transitionTxt, \ selfTxt, \ 
            txtOffset, \ selfTxtOffset, \ ftSize, \ ftFamily, \ ftSlant, 
            markovTxt}, 
α\  = \ 
            getOption[AngularOffset, \ 0, \ opts]; 
r\  = \ 
            getOption[StateSize, \ 1/4, \ opts]; 
r2\  = \ 
            getOption[SelfCircleSize, \ 1/5, \ opts]; 
          r3\  = \ getOption[SelfCircleOffset, \ 1/3, \ 
              opts]; 

lt\  = \ 
            getOption[LineThickness, \ 1/2, \ 
              opts]; 
h\  = \ 
            getOption[BendAngle, 20\ Degree, \ 
              opts]; 
arrowAngle1\  = \ 
            getOption[ArrowAngleTop, 10\ Degree, \ 
              opts]; 
arrowAngle2\  = \ 
            getOption[ArrowAngleBottom, 10\ Degree, \ opts]; 
          arrowLength\  = \ 
            getOption[ArrowLength,  .075, \ 
              opts]; 
stateColor\  = \ 
            getOption2[StateColor, GrayLevel[ .75], opts]; 
          hmmStyle\  = \ getOption3[MarkovStyle, \ Circular, \ opts]; 
          txtOffset\  = \ getOption[TextOffset, 1.1, \ opts]; 
          selfTxtOffset\  = \ getOption[SelfTextOffset, 0.1, \ opts]; 
          ftSize\  = \ getOption3[FontSize, \ 12, \ opts]; 
          ftFamily\  = \ getOption3[FontFamily, \ "Times", \ opts]; 
          ftSlant\  = \ getOption3[FontSlant, \ "Italic", \ opts]; 
          markovTxt\  = \ 
            getOption3[MarkovText, \ True, \ 
              opts]; 

at\  = \ 
            Transpose[a]; 
n\  = \ 
            Length[a]; 

nonZero\  = \ 
            Table[If[\((a[\([i, j]\)]\  != \ 0)\) && \((i != j)\), \ {i, 
                  j}], {i, 1, n}, {j, 1, n}]; 
          nonZero\  = \ Partition[Flatten[eliminateNulls[nonZero]], 2]; 
          nonZeroSelf\  = \ 
            Table[If[a[\([i, i]\)]\  != 0, \ {i, i}], \ {i, 1, n}]; 
          nonZeroSelf\  = \ 
            eliminateNulls[nonZeroSelf]; 
centers\  = 
            If[hmmStyle\  == \ Circular, 
              Table[{Cos[2  Pi\ \((i/n)\)\  + \ α], \ 
                    Sin[2  Pi\ \((i/n)\)\  + \ α]}\  // \ N, \ {i, 0, 
                  n - 1}], \(-1\), \(-1\)]; 
          If[centers\  == \ \(-1\), 
            centers\  = \ 
              If[hmmStyle\  == \ Sequential, \ 
                Table[{i, 0}, {i, 0, n - 1}], \(-1\), \(-1\)]]; 
          If[centers\  == \ \(-1\), \ 
            centers\  = \ 
              hmmStyle]; 

linesTmp\  = \ 
            Table[If[
                at[\([i, j]\)]\  > \ 0\  && \ 
                  i\  != \ j, \ {centers[\([i]\)], \ 
                  centers[\([j]\)]}], \ {i, \ 1, \ n}, \ {j, \ 1, n}]; 
          lines\  = \ \(Partition[#, 2] &\)\  /@ 
              Partition[Flatten[eliminateNulls[linesTmp]], 4]; 
          arcPts\  = \ \(thirdPoint[#, \ h] &\)\  /@ \ 
              lines; 
arcs\  = \ 
            circleParameters\  /@ \ arcPts; 

states\  = \ 
            Partition[Flatten[lines], 
              2]; 
states\  = \ \({#, \ r} &\)\  /@ \ 
              states; 

candidates\  = \(\((circleIntersections\  @@ \ #)\) &\)\  /@ \(({states, \ 
                    Join\  @@ \ 
                      Table[{arcs[\([i]\)], \ arcs[\([i]\)]}, \ {i, \ 
                          1\ Length[arcs]}]}\  // \ 
                  Transpose)\); 

third\  = \ 
            Join\  @@ \ Table[{Last[arcPts[\([i]\)]], \ 
                  Last[arcPts[\([i]\)]]}, {i, \ 1, \ 
                  Length[arcPts]}]; 
intersect\  = 
            Partition[
              minDistance\  /@ \ \(({candidates, \ third}\  // Transpose)\), 
              2]; 

arcsComplete\  = \((\(\((arcAngles\ \  @@ #)\) &\)\  /@ \ Transpose[{arcs, intersect}])\); 
          thetas\  = \ \(Last[Flatten[#]] &\)\  /@ \ 
              arcsComplete; 

arrows\  = \ \(Polygon[ArrowHead[{arrowAngle1, \ arrowAngle2, \ arrowLength}, \ 
                    convert[#]]] &\)\  /@ \ \(({arcs, \ thetas}\  // \ 
                  
                  Transpose)\); 

selfCenters\  = \ 
If[hmmStyle\  == \ Circular, 
              Table[{\((1\  + \ r3)\) 
                      Cos[2  Pi\ \((i/n)\)\  + \ α], \ \((1 + r3)\) 
                      Sin[2  Pi\ \((i/n)\)\  + \ α]}\  // \ N, \ {i, 
                  0, n - 1}], \(-1\), \(-1\)]; 
If[
            selfCenters\  == \ \(-1\), \ 
            selfCenters\  = \ 
If[
                hmmStyle\  == \ Sequential, \ 
                Table[{i + 1/1000, r3}\  // \ N, \ {i, 0, 
                    n - 1}], \(-1\), \(-1\)]]; 
			\ 
If[selfCenters\  == \ \(-1\), 
            selfCenters\  = \ \(\((# + {1/1000, r3})\) &\)\  /@ \ 
                centers]; 

selfArcs\  = \ Table[{selfCenters[\([i]\)], \ r2}, {i, 1, 
                n}]; 
		
intersect2\  = 
            Table[circleIntersections[{centers[\([i]\)], 
                  r}, \ {selfCenters[\([i]\)], \ r2}], {i, 1, 
                n}]; 
				
selfComplete\  = \((\(\((selfAngles\ \  @@ #)\) &\)\  /@ \ 
                Transpose[{selfArcs, intersect2}])\); 
          selfComplete\  = \ 
            Table[If[a[\([i, i]\)]\  != \ 0, selfComplete[\([i]\)]], {i, 1, 
                n}]; 
selfComplete\  = \ 
            eliminateNulls[
              selfComplete]; 

thetas2\  = \ \(Last[Flatten[#]] &\)\  /@ \ 
              selfComplete; 

selfArcs\  = \ Table[If[a[\([i, i]\)]\  != \ 0, selfArcs[\([i]\)]], {i, 1, n}]; 
          selfArcs\  = \ 
            eliminateNulls[
              selfArcs]; 

arrows2\  = \ \(Polygon[
                  ArrowHead[{arrowAngle1, \ arrowAngle2, \ arrowLength}, \ 
                    convert[#]]] &\)\  /@ \ \(({selfArcs, \ thetas2}\  // \ 
                  Transpose)\); 
	
stateTxt\  = 
            Table[Text[Subscript["S", ToString[i]], centers[\([i]\)], \ 
                Background -> GrayLevel[1. ]], {i, 1, \ n}]; 
          transitionTxt\  = \ \(TextConvert[#, 
                  txtOffset] &\)\  /@ \ \(({arcsComplete, nonZero}\  // \ 
                  Transpose)\); 
          selfTxt\  = \ \(TextConvert[#, 
                  selfTxtOffset] &\)\  /@ \(({selfComplete, 
                    nonZeroSelf}\  // \ 
                  Transpose)\); 
	
g0\  = \ 
            Graphics[\({stateColor, Disk[#, r]} &\)\  /@ \ centers]; 
          g1\  = Graphics[\({AbsoluteThickness[lt], 
                    Circle[#, \ r]} &\)\  /@ \ centers]; 
          g2\  = \ Graphics[{AbsoluteThickness[
                  lt], \(\((Circle\  @@ \ #)\) &\)\  /@ \ arcsComplete}]; 
          g3\  = \ Graphics[{AbsoluteThickness[ .1], \ arrows}]; 
          g4\  = \ Graphics[{AbsoluteThickness[
                  lt], \(\((Circle\  @@ \ #)\) &\)\  /@ \ selfComplete}]; 
          g5\  = \ Graphics[{AbsoluteThickness[ .1], \ arrows2}]; 
          g6\  = \ Graphics[{stateTxt, \ transitionTxt, 
                selfTxt}]; 

 (*\ 
            without\ arrows\ *) 

glist\  = \ If[markovTxt\  == \ True, \ {g4, g0, g1, \ g2, g3, g5, \ g6}, \ {g0, g1, 
                g2, g3, g4, g5}, {g0, g1, g2, g3, g4, g5, 
                g6}]; 
		
Show[glist, 
            TextStyle\  -> {FontSize -> ftSize, \ FontFamily -> ftFamily, \ 
                FontSlant -> ftSlant}, AspectRatio -> Automatic, \ 
            Frame -> False, \ PlotRange -> All, 
            ImageSize -> 400]
];\)\)

Hidden Markov model drawing function (DrawHMM)

Description

DrawHMM[A, B, opts] takes the state transition matrix   A, observable matrix B zero, and zero or more options opt, and draws the corresponding hidden Markov model. When the Markov model is labeled, some subscripts may be reversed. In that case,  FixText[Graphic, {{i1,j1},{i2,j2},...}] will switch the {ik,jk} subscripts. The function SubstituteValues[Graphic, A] will substitute actual values of the matrix A into the Markov model plot.

Options:

* HMMStyle -> {Circular, Sequential}: arrange states  in a circle or in a line.
* AngularOffset -> 0 Degree:  offset for circularly arranged Markov models.
* StateSize -> 1/4: relative size of state circles.
* StateColor -> GrayLevel[.75]: color of each state.
* LineThickness -> 1/2 (point): default line thickness  for drawing.
* BendAngle -> 20 Degree:  curvature of state-transition arcs.
* SelfCircleSize -> 1/5: relative radius of self-connections.
* SelfCircleOffset -> 1/3: relative offset of self-connections from center of states.
* ArrowAngleTop -> 10 Degree: arrow tip spread angle.
* ArrowAngleBottom -> 10 Degree: arrow tail  angle.
* ArrowLength -> 0.075: relative length of arrows.
* TextOffset -> 1.1: transition arcs offset for text labeling.
* SelfTextOffset -> 0.1: self-transition arcs offset for text labeling.
* MarkovText -> True: whether to label Markov model with text.
* FontSize -> 12: what font size to use.
* FontFamily -> "Times": what font family to use.
* FontSlant -> "Italic": what font slant to use.
* ObservableStyle -> {Style1, Style2,...}: 
* StateText -> True: whether to label states.

Definition

\(\(DrawHMM[a_, \ b_, opts___] := 
        Module[{at, \ bt, \ c, \ n, \ α, \ r, \ lt, \ g0, g1, g2, g3, 
            g4, g5, g6, \ glist, lines, linesTmp, \ h, \ arcPts, \ arcs, \ 
            states, \ candidates, third, \ intersect, arcsComplete, \ 
            arrows, \ arrowAngle1, \ arrowAngle2, arrowLength, \ 
            selfCenters, \ r2, \ r3, \ intersect2, \ selfArcs, 
            selfComplete, \ thetas2, \ arrows2, \ stateColor, \ hmmStyle, \ 
            stateTxt, nonZero, \ nonZeroSelf, \ transitionTxt, \ selfTxt, \ 
            txtOffset, \ selfTxtOffset, \ ftSize, \ ftFamily, \ ftSlant, 
            markovTxt, obs, \ obsStyle, obsVec, dskList, \ stateTxtFlag, 
            txtList}, α\  = \ 
            getOption[AngularOffset, \ 0, \ opts]; 
r\  = \ 
            getOption[StateSize, \ 1/4, \ opts]; 
r2\  = \ 
            getOption[SelfCircleSize, \ 1/5, \ 
              opts]; 
r3\  = \ 
            getOption[SelfCircleOffset, \ 1/3, \ 
              opts]; 
	
lt\  = \ 
            getOption[LineThickness, \ 1/2, \ 
              opts]; 
h\  = \ 
            getOption[BendAngle, 20\ Degree, \ opts]; 
          arrowAngle1\  = \ getOption[ArrowAngleTop, 10\ Degree, \ opts]; 
          arrowAngle2\  = \ getOption[ArrowAngleBottom, 10\ Degree, \ opts]; 
          arrowLength\  = \ getOption[ArrowLength,  .075, \ opts]; 
          stateColor\  = \ getOption2[StateColor, GrayLevel[ .75], opts]; 
          hmmStyle\  = \ getOption3[HMMStyle, \ Circular, \ opts]; 
          txtOffset\  = \ getOption[TextOffset, 1.1, \ opts]; 
          selfTxtOffset\  = \ getOption[SelfTextOffset, 0.1, \ opts]; 
          ftSize\  = \ getOption3[FontSize, \ 12, \ opts]; 
          ftFamily\  = \ getOption3[FontFamily, \ "Times", \ opts]; 
          ftSlant\  = \ getOption3[FontSlant, \ "Italic", \ opts]; 
          markovTxt\  = \ getOption3[MarkovText, \ True, \ opts]; 
          obsStyle\  = \ 
            getOption3[
              ObservableStyle, \ {RGBColor[1, 0, 0], \ RGBColor[0, 1, 0], \ 
                RGBColor[0, 0, 1], \ RGBColor[1, 1, 0], \ 
                RGBColor[1, 0, 1], \ RGBColor[0, 1, 1]}, \ opts]; 
          stateTxtFlag\  = \ 
            getOption3[StateText, \ True, \ 
              opts]; 
		
at\  = \ 
            Transpose[a]; 
bt\  = \ 
            Transpose[b]; 
n\  = \ 
            Length[a]; 
obs\  = \ 
            Length[b]; 
obsStyle\  = \ 
            Join\  @@ Table[
                obsStyle, \ {Ceiling[
                    obs/Length[
                        obsStyle]]}]; 
	
nonZero\  = \ 
            Table[If[\((a[\([i, j]\)]\  != \ 0)\) && \((i != j)\), \ {i, 
                  j}], {i, 1, n}, {j, 1, n}]; 
          nonZero\  = \ 
            Partition[Flatten[eliminateNulls[nonZero]], 
              2]; 
		
nonZeroSelf\  = \ 
            Table[If[a[\([i, i]\)]\  != 0, \ {i, i}], \ {i, 1, n}]; 
          nonZeroSelf\  = \ 
            eliminateNulls[nonZeroSelf]; 
		
centers\  = 
            If[hmmStyle\  == \ Circular, 
              Table[{Cos[2  Pi\ \((i/n)\)\  + \ α], \ 
                    Sin[2  Pi\ \((i/n)\)\  + \ α]}\  // \ N, \ {i, 0, 
                  n - 1}], Table[{i, 0}, {i, 0, n - 1}], 
              Table[{i, 0}, {i, 0, n - 1}]]; 
          linesTmp\  = \ 
            Table[If[
                at[\([i, j]\)]\  > \ 0\  && \ 
                  i\  != \ j, \ {centers[\([i]\)], \ 
                  centers[\([j]\)]}], \ {i, \ 1, \ n}, \ {j, \ 1, n}]; 
          lines\  = \ \(Partition[#, 2] &\)\  /@ 
              Partition[Flatten[eliminateNulls[linesTmp]], 4]; 
          arcPts\  = \ \(thirdPoint[#, \ h] &\)\  /@ \ 
              lines; 
arcs\  = \ 
            circleParameters\  /@ \ 
              arcPts; 

states\  = \ 
            Partition[Flatten[lines], 
              2]; 
states\  = \ \({#, \ r} &\)\  /@ \ 
              states; 
		
candidates\  = \(\((circleIntersections\  @@ \ #)\) &\)\  /@ \(({states, \ 
                    Join\  @@ \ 
                      Table[{arcs[\([i]\)], \ arcs[\([i]\)]}, \ {i, \ 
                          1\ Length[arcs]}]}\  // \ 
                  Transpose)\); 

third\  = \ 
            Join\  @@ \ Table[{Last[arcPts[\([i]\)]], \ 
                  Last[arcPts[\([i]\)]]}, {i, \ 1, \ Length[arcPts]}]; 
          intersect\  = 
            Partition[
              minDistance\  /@ \ \(({candidates, \ third}\  // Transpose)\), 
              2]; arcsComplete\  = \((\(\((arcAngles\ \  @@ #)\) &\)\  /@ \ 
                Transpose[{arcs, intersect}])\); 
          thetas\  = \ \(Last[Flatten[#]] &\)\  /@ \ arcsComplete; 
          arrows\  = \ \(Polygon[
                  ArrowHead[{arrowAngle1, \ arrowAngle2, \ arrowLength}, \ 
                    convert[#]]] &\)\  /@ \ \(({arcs, \ thetas}\  // \ 
                  Transpose)\); 
          selfCenters\  = \ 
            If[hmmStyle\  == \ Circular, 
              Table[{\((1\  + \ r3)\) 
                      Cos[2  Pi\ \((i/n)\)\  + \ α], \ \((1 + r3)\) 
                      Sin[2  Pi\ \((i/n)\)\  + \ α]}\  // \ N, \ {i, 
                  0, n - 1}], \ 
              Table[{i + 1/1000, r3}\  // \ N, \ {i, 0, n - 1}], \ 
              Table[{i + 1/1000, \ r3}\  // \ N, \ {i, \ 0, \ 
                  n - 1}]]; 
selfArcs\  = \ 
            Table[{selfCenters[\([i]\)], \ r2}, {i, 1, 
                n}]; 
		
intersect2\  = 
            Table[circleIntersections[{centers[\([i]\)], 
                  r}, \ {selfCenters[\([i]\)], \ r2}], {i, 1, 
                n}]; 
				
selfComplete\  = \((\(\((selfAngles\ \  @@ #)\) &\)\  /@ \ 
                Transpose[{selfArcs, intersect2}])\); 
          selfComplete\  = \ 
            Table[If[a[\([i, i]\)]\  != \ 0, selfComplete[\([i]\)]], {i, 1, 
                n}]; selfComplete\  = \ 
            eliminateNulls[
              selfComplete]; 

thetas2\  = \ \(Last[Flatten[#]] &\)\  /@ \ selfComplete; 
          selfArcs\  = \ 
            Table[If[a[\([i, i]\)]\  != \ 0, selfArcs[\([i]\)]], {i, 1, n}]; 
          selfArcs\  = \ 
            eliminateNulls[
              selfArcs]; 

arrows2\  = \ \(Polygon[ArrowHead[{arrowAngle1, \ arrowAngle2, \ arrowLength}, \ 
                    convert[#]]] &\)\  /@ \ \(({selfArcs, \ thetas2}\  // \ 
                  Transpose)\); 

stateTxt\  = \ 
            Table[Text[Subscript["S", ToString[i]], centers[\([i]\)], \ 
                Background -> GrayLevel[1. ]], {i, 1, \ 
                n}]; 
		transitionTxt\  = \ \(TextConvert[#, 
                  txtOffset] &\)\  /@ \ \(({arcsComplete, nonZero}\  // \ 
                  Transpose)\); 
		selfTxt\  = \ \(TextConvert[#, 
                  selfTxtOffset] &\)\  /@ \ \(({selfComplete, 
                    nonZeroSelf}\  // \ 
                  Transpose)\); 
		
obsVec\  = \ 
            TransformVector /@ \ bt; 
          dskList\  = \ 
            Table[DiskList[centers[\([i]\)], \ r, \ obsVec[\([i]\)], \ 
                obsStyle], {i, 1, n}]; 
          txtList\  = 
            If[stateTxtFlag\  == \ True, {stateTxt, \ transitionTxt, 
                selfTxt}, \ {transitionTxt, selfTxt}\ , {transitionTxt, 
                selfTxt}]; 

g0\  = \ 
            Graphics[dskList]; 
          g1\  = Graphics[\({AbsoluteThickness[lt], 
                    Circle[#, \ r]} &\)\  /@ \ centers]; 
          g2\  = \ Graphics[{AbsoluteThickness[
                  lt], \(\((Circle\  @@ \ #)\) &\)\  /@ \ arcsComplete}]; 
          g3\  = \ Graphics[{AbsoluteThickness[ .1], \ 
                arrows}]; 
g4\  = \ 
            Graphics[{AbsoluteThickness[
                  lt], \(\((Circle\  @@ \ #)\) &\)\  /@ \ selfComplete}]; 
          g5\  = \ Graphics[{AbsoluteThickness[ .1], \ 
                arrows2}]; 
g6\  = \ 
            Graphics[txtList]; 

glist\  = \ 
            If[markovTxt\  == \ True, \ {g0, g1, \ g2, g3, g4, g5, 
                g6}, \ {g0, g1, g2, g3, g4, g5}, {g0, g1, g2, g3, g4, 
                g5}]; 

 (*\ \(glist\  = \ If[markovTxt\  == \ True, \ {g0, g1, g2, g4, g5}, \ {g0, g1, g2, g4}, {g0, 
                    g1, g2, 
                    g4}];\)\ *) 

Show[
            glist, TextStyle\  -> {FontSize -> ftSize, \ 
                FontFamily -> ftFamily, \ FontSlant -> ftSlant}, 
            AspectRatio -> Automatic, \ Frame -> False, 
            PlotRange -> All, \ 
            ImageSize -> 400]
];\)\)

Draw observation sequence (DrawObservationSequence)

DrawObservationSequence[o_,  styles_, 
          opts___] := 
Module[{g1, g2,  length}, 
          length = 40; 
g1 = 
            Table[{styles[[o[[i]]]], 
                Disk[{Mod[i - 1, 
                      length], -Quotient[i - 1, length]},  .4]}, {i, 1, 
                Length[o]}]; 
g2 = 
            Table[Circle[{Mod[i - 1, 
                    length], -Quotient[i - 1, length]},  .4], {i, 1, 
                Length[o]}]; 
Show[Graphics[{g1, g2}], 
            AspectRatio -> Automatic, ImageSize -> 400, 
            opts]
];

Draw observation sequence transpose (DrawObservationSequenceTranspose)

DrawObservationSequenceTranspose[o_,  styles_, 
          opts___] := 
Module[{g1, g2,  length}, 
          length = 40; 
g1 = 
            Table[{styles[[o[[i]]]], 
                Disk[{-Quotient[i - 1, length], 
                    Mod[i - 1, length]},  .4]}, {i, 1, 
                Length[o]}]; 
g2 = 
            Table[Circle[{-Quotient[i - 1, length], 
                  Mod[i - 1, length]},  .4], {i, 1, 
                Length[o]}]; 
Show[Graphics[{g1, g2}], 
            AspectRatio -> Automatic, ImageSize -> 400, 
            opts]
];

Import definitions

Read HMM file (ReadHMM)

ReadHMM[
          file_] := 
Module[{A, B, π,  tmp,  
            nStates,  nObservables}, 
fp  =  
            OpenRead[file]; 
tmp  =  
            Read[fp,  Word]; 
nStates  =  
            Read[fp,  Number]; 
tmp = 
            Read[fp, Word]; 
nObservables = 
            Read[fp, Number]; 
tmp = 
            Read[fp, Word]; 
A = 
            Read[fp, 
              Table[Number, {nStates}, {nStates}]]; 
tmp = 
            Read[fp, Word]; 
B = 
            Read[fp, 
              Table[Number, {nObservables}, {nStates}]]; 
tmp = Read[fp, Word]; 
π = 
            Read[fp, Table[Number, {nStates}]]; 
Close[
            fp]; 
{A, 
            B, π}
];

Read observation sequence file (ReadObservationSequences)

ReadSingleObservationSequence[
          fp_] := 
Module[{o, T, 
            label}, 
T = 
            Read[fp, Number]; 
label = 
            Read[fp, Word]; 
o = 
            Read[fp, Table[Number, {T}]]; 
o + 
            1
];

ReadObservationSequences[
          file_] := 
Module[{fp, M, 
            os}, 
fp  =  
            OpenRead[file]; 
M = 
            Read[fp, Number]; 
Read[fp, 
            Number]; 
os = 
            Table[ReadSingleObservationSequence[
                fp], {M}]; 
Close[
            fp]; 
os
];

HMM definitions

Floor HMM values (NormalizeHMM, FloorHMM)

NormalizeRows[
          mat_] := 
(#/(Plus  @@  #)) &  /@  mat;

NormalizeColumns[mat_] := 
Transpose[
          NormalizeRows[mat  //  Transpose]];

NormalizeHMM[λ_] := 
Module[{a, 
            b, π}, 
{a, 
              b, π} = λ; 
{NormalizeRows[a], 
              NormalizeColumns[b], NormalizeRows[{π}] // First}  //  
            N];

FloorHMM[λ_,  which_: {1, 1, 1},  
          floor_:  10^-9] := 
Module[{a, b, π,  
            af,  bf,  πf, n, L,  
            individualFloors}, 
individualFloors = 
            which*floor; 
{a, 
              b, π} = λ; 
n = 
            Length[a]; 
L = 
            Length[b]; 
af = 
            Table[Max[a[[i, j]], individualFloors[[1]]], {i, 1, 
                n}, {j, 1, n}]; 
bf = 
            Table[Max[b[[k, j]], individualFloors[[2]]], {k, 1, 
                L}, {j, 1, n}]; 
πf = 
            Table[Max[π[[i]], individualFloors[[3]]], {i, 1, 
                n}]; 
NormalizeHMM[{af, 
              bf, πf}]
];

Random HMM model (RandomHMM)

RandomHMM[n_, L_, 
          structure_] := 
Module[{a, 
            b, π}, 
a = 
            If[structure == 1, 
              Table[If[j ≥ i,  Random[Real, {0, 1}], 0], {i, 
                  1, n}, {j, 1, n}], 
Table[
                Random[Real, {0, 1}], {n}, {n}], 
              Table[Random[Real, {0, 1}], {n}, {n}]]; 
b = 
            Table[Random[Real, {0, 1}], {L}, {n}]; 
π = 
            Table[Random[
                Real, {0, 1}], {n}]; 
NormalizeHMM[{a, 
              b, π}]
];

Generate Markov (and HMM) sequences (GenerateObservationSequence)

rr := Random[Real,  {0. ,  1. }];

TransformProbabilityVector[
          p_] := 
Table[
            Plus  @@  Take[p, {1, i}], {i, 1, Length[p]}]  //  
          N;

PickOne[ p_] := 
Module[{r}, 
          r  =  rr; (Position[(r  <=  #) &  /@  
                    TransformProbabilityVector[p],  True]  //  
                Flatten)  //  First];

GenerateNextState[q_,  A_] := 
PickOne[
          A[[q]]];

GenerateStateSequence[λ_,  
          n_] := 
Module[{A, 
            B, π}, 
{A,  
              B, π}  =  λ; 
NestList[
            GenerateNextState[#, A] &, PickOne[π], 
            n - 1]
];

GenerateOneObservation[q_, 
          B_] := 
PickOne[(B  //  
                Transpose)[[q]]];

GenerateObservationSequence[λ_,  
          n_] := 
Module[{A, B, π, 
            q}, 
{A,  
              B, π}  =  λ; 
q = 
            GenerateStateSequence[λ, 
              n]; 
GenerateOneObservation[#, 
                B] &  /@  q];

GenerateBothSequences[λ_,  
          n_] := 
Module[{A, B, π, 
            q}, 
{A,  
              B, π}  =  λ; 
q = 
            GenerateStateSequence[λ, 
              n]; 
o = GenerateOneObservation[#, 
                  B] &  /@  q; 
{q, o}];

Forward-backward algorithm (ForwardBackward)

ForwardBackward[
          o_, λ_] := 
Module[{a, b, π, n, T, 
            c, α,  αs,  β,  βs}, 
{a, b, π} = λ; 
n = 
            Length[a]; 
T = 
            Length[o]; 
Do[α[1, i] = π[[i]] 
                b[[o[[1]], i]], {i, 1, n}]; 
c[1] = 
            1. /(Sum[α[1, i], {i, 1, 
                    n}]); 
Do[αs[1, i] = 
              c[1] α[1, i], {i, 1, 
              n}]; 
Do[
Do[α[
                  t + 1, j] = 
                Sum[αs[t, i] a[[i, j]], {i, 1, n}] 
                  b[[o[[t + 1]], j]], {j, 1, 
                n}]; 
c[t + 1] = 
              1. /Sum[α[t + 1, i], {i, 1, 
                    n}]; 
Do[αs[t + 1, j] = 
                c[t + 1] α[t + 1, j], {j, 1, n}], {t, 1, 
              T - 1}]; 

Do[β[T, 
                i] = 1. , {i, 1, n}]; 
Do[βs[T, i] = 
              c[T] β[T, i], {i, 1, 
              n}]; 
Do[
Do[β[t, 
                  i] = Sum[
                  a[[i, j]] b[[o[[t + 1]], j]] βs[t + 1, 
                      j], {j, 1, n}], {i, 1, 
                n}]; 
Do[βs[t, i] = 
                c[t] β[t, i], {i, 1, n}], {t, T - 1, 
              1, -1}]; 
{Table[αs[t, i], {i, 1, 
                n}, {t, 1, T}], Table[βs[t, i], {i, 1, n}, {t, 1, T}], 
            Table[c[t], {t, 1, T}]}
];

Evaluation problem (LogProb, NormProb)

LogProb[
          o_, λ_] := 
Module[{αs, βs, 
            c, T}, 
T = 
            Length[o]; 
{αs, βs, c} = 
            ForwardBackward[o, λ]; 
-Sum[
              Log[10, c[[t]]], {t, 1, 
                T}]
];

NormProb[o_, λ_] := 
        10^(LogProb[o, λ]/Length[o]);

Decoding problem (DecodeHMM)

ArgMax[v_] := Position[v, Max[v]]  //  Flatten  //  
          First;

DecodeHMM[
          o_, λ_] := 
Module[{a, 
            b, π, ϕ, ψ, q, n, 
            T}, 
{a, 
              b, π} = λ; 
n = 
            Length[a]; 
T = 
            Length[o]; 
Do[ϕ[1, i] = 
              Log[π[[i]]]  +  Log[b[[o[[1]], i]]], {i, 1, 
              n}]; 
Do[ϕ[t, j] = 
              Max[Table[ϕ[t - 1, i]  +  Log[a[[i, j]]], {i, 1, 
                      n}]] + Log[b[[o[[t]], j]]], {t, 2, T}, {j, 1, 
              n}]; 
Do[ψ[t, j] = 
              ArgMax[Table[ϕ[t - 1, i] + Log[a[[i, j]]], {i, 1, 
                    n}]], {t, 2, T}, {j, 1, n}]; 
q[T] = 
            ArgMax[Table[ϕ[T, i], {i, 1, n}]]; 
Do[
            q[t] = ψ[t + 1, q[t + 1]], {t, T - 1, 
              1, -1}]; 
Table[
            q[t], {t, 1, T}]
];

Training problem (TrainHMM)

convHMM[param1_,  
          param2_] := 
Module[{df}, 
df  =  Flatten[param1]  -  
              Flatten[param2]; 
Sqrt[(Plus  @@  (df . 
                      df))/Length[df]]
];

TrainHMMOneStep[
          o_, λ_] := 
Module[{fb, α, β, c, a, b, π,  L,  T,  M, anew, bnew, 
            bselect, πnew}, 
{a, 
              b, π} = λ; 
n = 
            Length[a]; 
L = 
            Length[b]; 
M = 
            Length[o]; 
T = 
            Length  /@  
              o; 
fb = ForwardBackward[#, λ] &  /@   o; 
den = 
            Plus  @@ Table[
                Table[{α, β, c} = fb[[m]]; 
                  Sum[α[[i, t]] β[[i, t]]/
                        c[[t]], {t, 1, T[[m]] - 1}], {i, 1, n}], {m, 
                  1, M}]; 
anew = 
            Plus  @@ Table[
                Table[{α, β, c} = fb[[m]]; 
                  Sum[α[[i, t]] a[[i, 
                            j]] b[[o[[m, t + 1]], 
                            j]] β[[j, t + 1]], {t, 1, 
                        T[[m]] - 1}]/den[[i]], {i, 1, n}, {j, 1, 
                    n}], {m, 1, 
                  M}]; 

bselect = 
            Table[Table[
                Position[o[[m]], k]  //  Flatten, {k, 1, L}], {m, 1, 
                M}]; 
den = 
            den + (Plus @@ 
                  Table[Table[{α, β, c} = 
                        fb[[m]]; α[[i, 
                            T[[m]]]] β[[i, T[[m]]]]/
                          c[[T[[m]]]], {i, 1, n}], {m, 1, 
                      M}]); 
bnew = 
            Plus @@ Table[{α, β, c} = fb[[m]]; 
                Table[(Plus @@ ((α[[j, #]] β[[j, #]]/c[[#]]) & /@  bselect[[m, k]]))/den[[j]], {k, 1, 
                    L}, {j, 1, n}], {m, 1, M}]; 
If[
            M > 1, πnew = {1. }~Join~
                Table[0. , {n - 1}], {α, β, c} = 
              fb[[1]]; 
πnew = 
              Table[α[[i, 1]] β[[i, 1]]/c[[1]], {i, 
                  1, n}]]; 
NormalizeHMM[{anew, 
              bnew,  πnew}]
];

TrainHMM[o_, λinitial_, max_, 
          tst_:  .001] := 
NestWhileList[
          TrainHMMOneStep[
              o, #] &,  λinitial,  (convHMM[#1, #2] > tst) &,  
          2, max];