/*
 * Example from section 2 of "Time Recursive EM Examples",
 * R. Perry, June 1998
 */

/* load b4 and vit first */


K = 10;	/* sliding window size */


double em_r( double _yK)
{
/*
 * time-recursive EM update for y = c*b + z problem
 * using a sliding window.
 *
 * Given new value for yK:
 *  E-step: estimate new value for c (called d here)
 *  M-step: compute ML estimate of b_buf using viterbi algorithm
 */

  vector j;
  int i, jn1;
  double s, b1, y1;

  y_buf(K) = _yK;
  b_buf(K) = 1;

  for( i = 0; i < 3; ++i) {	/* fixed at 3 iterations */

    /* E-step */

	j = find( b_buf == 1);
	jn1 = cols(j);
	d = c;
	if( jn1 > 0) {
	  s = sum( y_buf(j));
	  d += (s - jn1 * c) / (n1 + jn1);
	}

    /* M-step */

        b_buf = vit( y_buf, 2, istate) - 1;

    if( verbose > 1)
      printf( "em_r: yK = %g, d = %g, jn1 = %g\n", _yK, d, jn1);
  }

  b1 = b_buf(1);
  y1 = y_buf(1);

  if( K > 1) {
    b_buf = [b_buf(2:K), 0];	/* shift */
    y_buf = [y_buf(2:K), 0];
  }

  if( b1 == 1) {
    ++n1;
    c += (y1 - c) / n1;
  }

  istate = b1 + 1;	/* initial viterbi state for next call */

  return b1;
}



void iter_r( void)
{
  int i, k;
  vector j;

  if( K > N) {
    print( "iter_r: K > N not allowed!");
    return;
  }

  c = 0;
  n1 = 0;
  istate = 0;		/* initial viterbi state unknown */

  y_buf = y(1:K);	/* initialize pipeline */
  b_buf = ones(1,K);

  for( i = K, k = 1; i <= N; ++i, ++k) {
      if( verbose > 1) printf( "k = %g\n", k);
      b(k) = em_r( y(i));
  }

  if( K > 1)
    b(k:N) = b_buf(1:K-1);	/* drain the pipeline */
}
