/*
 * Example from sections 1-4 of "Dual Binomial EM Examples",
 * R. Perry, June 1998
 */

s = 0.4;
n1 = 9;
n0 = 8;

q1 = 0.2;
q2 = 0.8;

void iter( int _n)
{
/*
 * perform n iterations of EM for dual binomial problem
 *
 * Given values for s, q1, and q2:
 * E-step: estimate values for n11, n12, n01, n02
 * M-step: compute ML estimates of q1 and q2
 */

  int i, n;
  double t, t1, t2;

  n = (narg > 0) ? _n : 1;

  for( i = 0; i < n; ++i) {

    /* E-step */

      if( q1 == 0 && q2 == 0)	/* degenerate case */

	n11 = n12 = 0;

      else {

	t1 = s*q1;
	t2 = (1-s)*q2;
	t  = t1 + t2;

	n11 = t1 * n1 / t;

	n12 = t2 * n1 / t;
      }

      if( q1 == 1 && q2 == 1)	/* degenerate case */

	n01 = n02 = 0;

      else {

        t1 = s*(1-q1);
        t2 = (1-s)*(1-q2);
        t  = t1 + t2;

        n01 = t1 * n0 / t;

        n02 = t2 * n0 / t;
      }

    /* M-step */

      q1 = n11/(n11+n01);

      q2 = n12/(n12+n02);

      p = s*q1 + (1-s)*q2;

      print(p);
  }

  print(q1,q2,n11,n12,n01,n02);
}

double p_ml( void)
{
/*
 * directly compute ML estimate of p
 */
  return n1/(n1+n0);
}
