2026.sagews

This is the Numerical Recipes 3rd edition code for an iterative improvement Hidden Markov Modeller naively ported to Mathematica.

(* data for the Hidden Markov Model should be stored in dtest *)

n = 50; (* number of levels of iterative improvement *)
nruns = 1; (* number of runs from initial random guesses *)
mstat = 3; (* number of states *)
kstates = 6; (* number of outcomes *)
BIG = 10^20; (* need this to handle the extremely low probabilities encountered *)
BIGI = N[1/BIG];
SeedRandom[1]; (* seed random numbers for reproducibility *)

tend = Dimensions[dtest][[1]];
jlindex = Table[Table[0, {i, 1, n}], {i, 1, nruns}];
For[run = 1, run <= nruns, run++,
 \[Alpha] = Table[Table[0, {i, 1, mstat}], {i, 1, tend}];
 \[Beta] = Table[Table[0, {i, 1, mstat}], {i, 1, tend}];
 pstate = \[Alpha];
 \[Alpha]norm = Table[0, {i, 1, tend}];
 \[Beta]norm = \[Alpha]norm;
 A = Table[Table[RandomReal[], {i, 1, mstat}], {i, 1, mstat}];
 b = Table[Table[RandomReal[], {i, 1, kstates}], {i, 1, mstat}];
 For[i = 1, i <= mstat, i++,
  A[[i, 1 ;; mstat]] = A[[i, 1 ;; mstat]]/Total[A[[i, 1 ;; mstat]]];
  b[[i, 1 ;; kstates]] = b[[i, 1 ;; kstates]]/
   Total[b[[i, 1 ;; kstates]]];
  ];
 (* Run the Hidden Markov Model *)
 powtab = Table[BIGI^(i - 6), {i, 0, 9}];
 For[irun = 1, irun <= n, irun++,
  \[Alpha] = 
   Table[Table[b[[ii, dtest[[i]]]], {ii, 1, mstat}], {i, 1, tend}];
  \[Beta] = Table[Table[1, {i, 1, mstat}], {i, 1, tend}];
  (* calculate \[Alpha] and \[Beta] and renormalize them - 
  from Numerical Recipes*)
  For[t = 2, t <= tend, t++,
   asum = 0;
   For[j = 1, j <= mstat, j++,
    sum = 0;
    For[i = 1, i <= mstat, i++,
     sum += \[Alpha][[t - 1, i]]*A[[i, j]]*b[[j, dtest[[t]]]];
     ];
    \[Alpha][[t, j]] = sum;
    asum += sum;
    ];
   \[Alpha]norm[[t]] = \[Alpha]norm[[t - 1]];
   If[asum < BIGI,
    ++\[Alpha]norm[[t]];
    For[j = 1, j <= mstat, j++,
     \[Alpha][[t, j]] *= BIG;
     ];
    ];
   ];
  \[Beta]norm[[tend]] = 0;
  For[t = tend - 1, t >= 1, t--,
   bsum = 0;
   For[i = 1, i <= mstat, i++,
    sum = 0;
    For[j = 1, j <= mstat, j++,
     sum += A[[i, j]]*b[[j, dtest[[t + 1]]]]*\[Beta][[t + 1, j]];
     ];
    \[Beta][[t, i]] = sum;
    bsum += sum;
    ];
   \[Beta]norm[[t]] = \[Beta]norm[[t + 1]];
   If[bsum < BIGI,
    ++\[Beta]norm[[t]];
    For[j = 1, j <= mstat, j++,
     \[Beta][[t, j]] *= BIG;
     ];
    ];
   ];
  lhood = 0;
  For[i = 1, i <= mstat, i++,
   lhood += \[Alpha][[1, i]]*\[Beta][[1, i]];
   ];
  lnorm = \[Alpha]norm[[1]] + \[Beta]norm[[1]];
  While[lhood < BIGI,
   lhood *= BIG;
   lnorm++;
   ];
  jlindex[[run, irun]] = Log[lhood] + lnorm*Log[BIGI];
  For[t = 1, t <= tend, t++,
   sum = 0;
   For[i = 1, i <= mstat, i++,
    pstate[[t, i]] = \[Alpha][[t, i]]*\[Beta][[t, i]];
    sum += pstate[[t, i]];
    ];
   (* sum=lhood*BIGI^(lnorm-\[Alpha]norm[[t]]-\[Beta]norm[[t]]); *)
   For[i = 1, i <= mstat, i++,
    pstate[[t, i]] *= 1/sum;
    ];
   ];
  bnew = Table[Table[0, {ii, 1, kstates}], {i, 1, mstat}];
  For[i = 1, i <= mstat, i++,
   denom = 0;
   For[t = 1, t <= tend, t++,
    term = (\[Alpha][[t, i]]*\[Beta][[t, i]])/lhood*
      powtab[[\[Alpha]norm[[t]] + \[Beta]norm[[t]] - lnorm + 6]];
    (*term=(\[Alpha][[t,i]]*\[Beta][[t,i]])/lhood*BIGI^(\[Alpha]norm[[
    t]]+\[Beta]norm[[t]]-lnorm);*)
    denom += term;
    bnew[[i, dtest[[t]]]] += term;
    ];
   For[j = 1, j <= mstat, j++,
    num = 0;
    For[t = 1, t <= tend - 1, t++,
     num += \[Alpha][[t, i]]*
        b[[j, dtest[[t + 1]]]]*\[Beta][[t + 1, j]]*
        powtab[[\[Alpha]norm[[t]] + \[Beta]norm[[t + 1]] - lnorm + 
           6]]/lhood;
     (*num+=\[Alpha][[t,i]]*b[[j,dtest[[t+1]]]]*\[Beta][[t+1,j]]*
     BIGI^(\[Alpha]norm[[t]]+\[Beta]norm[[t+1]]-lnorm);*)
     ];
    num /= lhood;
    A[[i, j]] *= num/denom;
    ];
   sumhold = Total[A[[i, 1 ;; mstat]]];
   For[j = 1, j <= mstat, j++,
    A[[i, j]] /= sumhold;
    ];
   For[k = 1, k <= kstates, k++,
    bnew[[i, k]] *= 1/denom;
    ];
   ];
  b = bnew;
  ];
 ]

Traceback (most recent call last): File "<stdin>", line 1, in <module> File "_sage_input_4.py", line 135, in <module> ]''', '/sagenb/sagenb/sage_notebook.sagenb/home/sdemjanenko/3/cells/2') File "/usr/local/sage/local/lib/python2.6/site-packages/sagenb-0.7.5.3-py2.6.egg/sagenb/misc/support.py", line 470, in syseval system.chdir(dir) File "/usr/local/sage/local/lib/python2.6/site-packages/sage/interfaces/mathematica.py", line 473, in chdir self.eval('SetDirectory["%s"]'%dir) File "/usr/local/sage/local/lib/python2.6/site-packages/sage/interfaces/mathematica.py", line 392, in eval s = Expect.eval(self, code, **kwds) File "/usr/local/sage/local/lib/python2.6/site-packages/sage/interfaces/expect.py", line 983, in eval return '\n'.join([self._eval_line(L, **kwds) for L in code.split('\n') if L != '']) File "/usr/local/sage/local/lib/python2.6/site-packages/sage/interfaces/mathematica.py", line 440, in _eval_line allow_use_file=allow_use_file, wait_for_prompt=wait_for_prompt) File "/usr/local/sage/local/lib/python2.6/site-packages/sage/interfaces/expect.py", line 634, in _eval_line return self._eval_line_using_file(line) File "/usr/local/sage/local/lib/python2.6/site-packages/sage/interfaces/expect.py", line 621, in _eval_line_using_file s = self._eval_line(self._read_in_file_command(tmp_to_use), allow_use_file=False) File "/usr/local/sage/local/lib/python2.6/site-packages/sage/interfaces/mathematica.py", line 440, in _eval_line allow_use_file=allow_use_file, wait_for_prompt=wait_for_prompt) File "/usr/local/sage/local/lib/python2.6/site-packages/sage/interfaces/expect.py", line 637, in _eval_line self._start() File "/usr/local/sage/local/lib/python2.6/site-packages/sage/interfaces/expect.py", line 472, in _start raise RuntimeError, "Unable to start %s"%self.__name RuntimeError: Unable to start mathematica