0
$\begingroup$

I've also posted this question here, since I'm not sure which community is more adequate for this type of question.

I've programmed imperatively the Baum-Welch algorithm for a 2 state hidden Markov model (HMM), as is explained in this book. When I try it with a size sample of 100, I usually get a decent estimate after 18 seconds. However, when I use a sample size of 500, it takes me between 290-350 seconds to get an improved estimation of the HMM parameters.

I get the impression the algorithm should not perform this badly. However I'm not sure why it does.

One thing I've discovered is that Mathematica is more functional oriented, so when I program iteratively I'm losing some computation 'power'. But is this really enough for this bad performance??

What is your experience with this algorithm?

I'm using as stopping criterion when the norm of the difference between the current estimate and the previous one is less than $10^{-3}$. This criterion just seemed simple enough. I'm also trying to use other stopping criteria, but it doesn't seems to alter much. The problem is each iteration of the algorithm, in a sample size of 500, takes almost 20 seconds more or less... This worries me.

Any help would be appreciated.

Here is my code.

EM1[y_, ite_] := 
Module[{n, listainicial, teta0, teta1, teta2, likeham, likehamtotal,
 pimattotal, sigmatotal, tetatotal, θk1, θk2, 
listak, listbk, fikyj, iterations},
n = Length[y] - 1;
listainicial = N[Table[0, {j, 1, n}]];
pimattotal = listainicial;
tetatotal = listainicial;
sigmatotal = listainicial;
likehamtotal = listainicial;

likeham[pimat_, fiyj_] := Module[{pizero, Pmat, counter, mult},
 pizero = pimat[[1]];
 Pmat = {pimat[[2]], pimat[[3]]};
 counter = 1;
 mult = IdentityMatrix[2];
 While[counter <= n,
  mult = 
   Dot[mult, 
    Dot[Pmat, {{fiyj[[1, counter]], 0}, {0, fiyj[[2, counter]]}}]];
  counter = counter + 1;
  ];
 Dot[Dot[pizero, mult], {1, 1}]
 ];

iterations = 1;
While[iterations <= ite,
(*initialization*)

pim0 = RandomVariate[UniformDistribution[{0, 1}]];
pim1 = RandomVariate[UniformDistribution[{0, 1}]];
pim2 = RandomVariate[UniformDistribution[{0, 1}]];

pimatk = N[{{pim0, 1 - pim0}, {pim1, 1 - pim1}, {pim2, 1 - pim2}}];
fikyj = N[Table[Table[1/n, {j, 1, n}], {i, 1, 2}]];

ak = N[Table[Table[0.1, {j, 1, n}], {i, 1, 2}]];
bk = N[Table[Table[1, {j, 1, n}], {i, 1, 2}]];
θk1 = {0.1, 0.1, 0.1};
θk2 = {0.5, 0.5, 0.5};

listak = listainicial;
listbk = listainicial;


While[Norm[θk1 - θk2] > 10^-3,

 θk1 = θk2;

 ak[[1, 1]] = pimatk[[1, 1]]*fikyj[[1, 1]];
 ak[[2, 1]] = pimatk[[1, 2]]*fikyj[[2, 1]];

  ak[[1, 1]] = ak[[1, 1]];
  ak[[2, 1]] = ak[[2, 1]];

 (*Induction*)
 counterE = 1;
 While [counterE <= n - 1,
  ak[[1, 
    counterE + 
     1]] = (pimatk[[1 + 1, 1]]*ak[[1, counterE]] + 
      pimatk[[2 + 1, 1]]*ak[[2, counterE]])*
    fikyj[[1, counterE + 1]];
  ak[[2, 
    counterE + 
     1]] = (pimatk[[1 + 1, 2]]*ak[[1, counterE]] + 
      pimatk[[2 + 1, 2]]*ak[[2, counterE]])*
    fikyj[[2, counterE + 1]];

   ak[[1, counterE + 1]] = ak[[1, counterE + 1]];
   ak[[2, counterE + 1]] = ak[[2, counterE + 1]];

  bk[[1, n - counterE]] = 
   pimatk[[1 + 1, 1]]*fikyj[[1, n - counterE + 1]]*
     bk[[1, n - counterE + 1]] + 
    pimatk[[1 + 1, 2]]*fikyj[[2, n - counterE + 1]]*
     bk[[2, n - counterE + 1]];
  bk[[2, n - counterE]] = 
   pimatk[[2 + 1, 1]]*fikyj[[1, n - counterE + 1]]*
     bk[[1, n - counterE + 1]] + 
    pimatk[[2 + 1, 2]]*fikyj[[2, n - counterE + 1]]*
     bk[[2, n - counterE + 1]];

  bk[[1, n - counterE]] = bk[[1, n - counterE]];
   bk[[2, n - counterE]] = bk[[2, n - counterE]];

  counterE = counterE + 1;
  ];

 τhijk = 
  Table[Table[
    Table[(ak[[h, j]]*pimatk[[h + 1, i]]*fikyj[[i, j + 1]]*
        bk[[i, j + 1]])/(Sum[
        Sum[ak[[hh, j]]*pimatk[[hh + 1, ii]]*fikyj[[ii, j + 1]]*
          bk[[ii, j + 1]], {ii, 1, 2}], {hh, 1, 2}]), {j, 1, 
      n - 1}], {i, 1, 2}], {h, 1, 2}];
 τijk = 
  Table[Table[
    ak[[i, j]]*
     bk[[i, j]]/(Sum[
        ak[[h, j]]*10^listak[[j]]*bk[[h, j]], {h, 1, 2}]), {j, 1, 
     n}], {i, 1, 2}];

 (*M-step*)

 teta0 = (Sum[
     Sum[(Sum[
          y[[j + 1]]*τijk[[i, j]], {j, 1, n}])*τijk[[i, 
         k]]*y[[k]]/(Sum[τijk[[i, j]], {j, 1, n}]) - 
       y[[k + 1]]*y[[k]]*τijk[[i, k]], {i, 1, 2}], {k, 1, 
      n}])/(Sum[
     Sum[(Sum[
           y[[j]]*τijk[[i, j]], {j, 1, 
            n}]/(Sum[τijk[[i, j]], {j, 1, n}]) - 
         y[[k]])*τijk[[i, k]]*y[[k]], {i, 1, 2}], {k, 1, n}]);

 teta1 = 
  Sum[τijk[[1, j]]*(y[[j + 1]] - teta0*y[[j]]), {j, 1, n}]/
   Sum[τijk[[1, j]], {j, 1, n}];

 teta2 = 
  Sum[τijk[[2, j]]*(y[[j + 1]] - teta0*y[[j]]), {j, 1, n}]/
   Sum[τijk[[2, j]], {j, 1, n}];

 sigmaest = 
  Sqrt[Sum[((y[[j + 1]] - teta1 - teta0*y[[j]])^2*τijk[[1, 
         j]] + (y[[j + 1]] - teta2 - teta0*y[[j]])^2*τijk[[2,
          j]]), {j, 1, n}]/(n)];

 θk2 = {teta0, teta1, teta2};

 (*Q[θ_,σ_]:=(τijk[[1,1]]*Log[pimatk[[1,
 1]]]+τijk[[2,1]]*Log[pimatk[[1,2]]])+Sum[Sum[
 Sum[τhijk[[h,i,j]]*Log[pimatk[[h+1,i]]],{j,1,n-1}],{i,1,
 2}],{h,1,2}]+Sum[Sum[τijk[[i,j]]*Log[1/(Sqrt[2*
 Pi]*σ)]*(-(y[[j+1]]-θ[[i+1]]-θ[[1]]*y[[
 j]])^2/(2*σ^2)),{j,1,n}],{i,1,2}];*)

 (*updates on remaining parameters*)
 pimatk[[1, 1]] = τijk[[1, 1]];
 pimatk[[1, 2]] = τijk[[2, 1]];(*π_ 0i=τijk[[i,
 1]]*)

 pimatk[[2, 1]] = 
  Sum[τhijk[[1, 1, j]], {j, 1, n - 1}]/
   Sum[τijk[[1, j]], {j, 2, n}];
 pimatk[[2, 2]] = 
  Sum[τhijk[[1, 2, j]], {j, 1, n - 1}]/
   Sum[τijk[[1, j]], {j, 2, n}];
 pimatk[[3, 1]] = 
  Sum[τhijk[[2, 1, j]], {j, 1, n - 1}]/
   Sum[τijk[[2, j]], {j, 2, n}];
 pimatk[[3, 2]] = 
  Sum[τhijk[[2, 2, j]], {j, 1, n - 1}]/
   Sum[τijk[[2, j]], {j, 2, n}];(*π_hi=(\!\(
\*UnderoverscriptBox[\(∑\), \(j = 2\), \(n\)]
\*SubscriptBox[\(τ\), \(hij\)]\))/\!\(
\*UnderoverscriptBox[\(∑\), \(j = 2\), \(n\)]
\*SubscriptBox[\(τ\), \(hj\)]\)*)
 (*Print["pimatk =",pimatk];*)

 counterM = 1;
 While[counterM <= n,
  fikyj[[1, counterM]] = (1/(Sqrt[2*Pi]*sigmaest))*
    E^(-(y[[counterM + 1]] - θk2[[2]] - θk2[[1]]*
            y[[counterM]])^2/(2*sigmaest^2));

  fikyj[[2, counterM]] = (1/(Sqrt[2*Pi]*sigmaest))*
    E^(-(y[[counterM + 1]] - θk2[[3]] - θk2[[1]]*
            y[[counterM]])^2/(2*sigmaest^2));
  counterM = counterM + 1;

  ];


 ];

likehamtotal[[iterations]] = likeham[pimatk, fikyj];
pimattotal[[iterations]] = pimatk;
tetatotal[[iterations]] = θk2;
sigmatotal[[iterations]] = sigmaest;

iterations = iterations + 1;
];

iterations = 1;
pos = 1;
While[iterations <= ite ,
If[likehamtotal[[iterations]] > likehamtotal[[pos]],
 pos = iterations];
iterations = iterations + 1
];
{pimattotal[[pos]], tetatotal[[pos]], sigmatotal[[pos]]}
];
$\endgroup$
  • 2
    $\begingroup$ Can you post your code? We can't really help you otherwise! $\endgroup$ – dr.blochwave Oct 1 '14 at 9:14
  • $\begingroup$ @blochwave Thanks for your interest. I've just posted my code. $\endgroup$ – An old man in the sea. Oct 1 '14 at 10:47
  • 2
    $\begingroup$ That is a lot of code for anyone on this site to have time to read/critique/analyze it fully and most likely, no one is going to. However from a quick glance at the code, I can most certainly say that the slowness is due to the procedural programming style. All those Sum[Sum[Sum[Sum... terms in nested Whiles are bound to give terrible performance. Writing the algorithm in a functional manner would certainly speed things up a lot. Unfortunately, it is hard (for me at least) to recommend specific changes to make to your code. If I were you, I'd start afresh with functional (use Igor's link) $\endgroup$ – rm -rf Oct 1 '14 at 14:37
  • $\begingroup$ I was about to comment on the nested Sum[Sum[Sum[...s! $\endgroup$ – dr.blochwave Oct 1 '14 at 14:40
  • $\begingroup$ I would really like to see a Compiled version of this though. I think it could be pretty fast. I don't see anything that can't be compiled to C. $\endgroup$ – C. E. Oct 1 '14 at 15:14
2
$\begingroup$

A very clear description of a Mathematica implementation of BW is given by Robert J Frey here. He seems to not have performance issues.

$\endgroup$
  • 1
    $\begingroup$ I Thank you for our interest. I'm definitely going to take a look at Frey's solution. However, it seems to rely mainly on functional programming. So, this doesn't answer my question. This is more a comment than an answer to my question. $\endgroup$ – An old man in the sea. Oct 1 '14 at 11:24
  • 3
    $\begingroup$ While this link may answer the question, it is better to include the essential parts of the answer here and provide the link for reference. Link-only answers can become invalid if the linked page changes. $\endgroup$ – m_goldberg Oct 1 '14 at 11:43
  • $\begingroup$ @m_goldberg I don't necessarily disagree, however here the OP knows what the BW algorithm does, and the implementation is nontrivial, so it seems cumbersome to reproduce it here. $\endgroup$ – Igor Rivin Oct 1 '14 at 11:58
  • 2
    $\begingroup$ Perhaps this would be better as a comment to the question. $\endgroup$ – m_goldberg Oct 1 '14 at 12:03
  • $\begingroup$ @IgorRivin Please participate here: meta.mathematica.stackexchange.com/q/1395/193 $\endgroup$ – Dr. belisarius Oct 1 '14 at 12:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.