C语言简单实现EM算法

#include <stdio.h>
#include <math.h>

int main()
{
  int i, j;
  //H正面:0；T反面:1.
  int n[5][10] = {0, 1, 1, 1, 0, 0, 1, 0, 1, 0,
                  0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
                  0, 1, 0, 0, 0, 0, 0, 1, 0, 0,
                  0, 1, 0, 1, 1, 1, 0, 0, 1, 1,
                  1, 0, 0, 0, 1, 0, 0, 0, 1, 0};
  int p = 5, q = 10; //p==行数 q==列数

  double A = 0.6;
  double sa = 0; //每一行的A的概率

  double B = 0.5;
  double sb = 0; //每一行的B的概率

  double fa = 0;
  double fb = 0;

  int count1 = 0; //反面的个数
  int count0 = 0; //正面的个数

  int flag = 1;
  int temp = 0;

  double c[100] = {0}; //将每一次得出的A的概率都保存在数组c中，方便比较前后两次的数值
  double d[100] = {0}; //将每一次得出的B的概率都保存在数组d中，方便比较前后两次的数值

  c[0] = A; //将自己设置的概率存入c[0]
  d[0] = B;

  //用于累加每一行中硬币正面向上和反面向上的个数
  double sa0 = 0, sa1 = 0;
  double sb0 = 0, sb1 = 0;

  while (flag)
  {
    for (i = 0; i < p; i++)
    {
      for (j = 0; j < q; j++)
      {
        if (n[i][j] == 1)
        {
          count1++;
        }
        else
        {
          count0++;
        }
      }
      //分别计算每行的A和B的概率
      sa = pow(A, count0) * pow(1 - A, count1);
      sb = pow(B, count0) * pow(1 - B, count1);
      
      //此行属于A或B的概率
      sa = sa / (sa + sb);
      //printf("%.2lf ", sa);
      sb = 1 - sa;
      //printf("%.2lf\n", sb);
      
      //正0，反1.
      //用于累加每一行中A硬币正面向上和反面向上的个数
      sa1 += count1 * sa;
      sa0 += count0 * sa;
      //用于累加每一行中B硬币正面向上和反面向上的个数
      sb1 += count1 * sb;
      sb0 += count0 * sb;

      sa = 0; //将概率置为0，为了下次重新计算
      sb = 0;
      count1 = 0; //将个数置为0，为了下次重新统计
      count0 = 0;
    }

    temp++; //用于计数一共循环了几次，同时用于将新的概率存入数组c和d中

    //计算出新的A和B的概率
    fa = sa0 / (sa0 + sa1);
    fb = sb0 / (sb0 + sb1);

    c[temp] = fa;
    d[temp] = fb;

    printf("第%d次:fa = %.7lf fb = %0.7lf\n", temp, fa, fb);

    for (i = 0; i < temp; i++)
    {
      if (c[i] == c[i + 1] && d[i] == d[i + 1])
      {
        flag = 0;
        printf("最终次：A = %.7lf B = %.7lf", c[i+1], d[i+1]);
        break;
      }
      if(fabs(c[i+1]-c[i]) <= 0.000001 && fabs(d[i+1]-d[i]) <= 0.000001)
      {
        flag = 0;
        printf("最终次：A = %.7lf B = %.7lf", c[i], d[i]);
        break;
      }
    }
    A = fa;
    B = fb;
    fa = 0;
    fb = 0;
    sa1 = 0;
    sa0 = 0;
    sb1 = 0;
    sb0 = 0;
  }

  return 0;
}

运行结果：