序列比对（九）从掷骰子说起HMM

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

typedef char State;
typedef char Result;
State state[] = {'F', 'L'};   // 所有的可能状态
Result result[] = {'1', '2', '3', '4', '5', '6'};   // 所有的可能符号
double init[] = {0.9, 0.1};    // 初始状态的概率向量
double emission[][6] = {   // 发射矩阵：行对应着状态，列对应着符号
  1.0/6, 1.0/6, 1.0/6, 1.0/6, 1.0/6, 1.0/6,
  0.1, 0.1, 0.1, 0.1, 0.1, 0.5
};
double trans[][2] = {   // 转移矩阵：行和列都是状态
  0.95, 0.05,
  0.1, 0.9
};
const int nstate = 2;
const int nresult = 6;
State* rst;   // 一串随机状态序列
Result* rres;  // 一串随机符号序列

int random(double* prob, const int n);
void randSeq(State* st, Result* res, const int n);
void printState(State* st, const int n);
void printResult(Result* res, const int n);

int main(void) {
  int i;
  int n = 300;
  int ll = 60;
  int nl, nd;
  if ((rst = (State*) malloc(sizeof(State) * (n))) == NULL || \
      (rres = (Result*) malloc(sizeof(Result) * (n))) == NULL) {
    fputs("Error: out of space!\n", stderr);
    exit(1);
  }
  randSeq(rst, rres, n);
  nl = n / ll;
  nd = n % ll;
  for (i = 0; i < nl; i++) {
    printf("Rolls\t");
    printResult(rres + ll * i, ll);
    printf("Die\t");
    printState(rst + ll * i, ll);
    printf("\n");
  }
  if (nd > 0) {
    printf("Rolls\t");
    printResult(rres + ll * i, nd);
    printf("Die\t");
    printState(rst + ll * i, nd);
    printf("\n");
  }
  free(rst);
  free(rres);
}

// 根据一个概率向量从0到n-1随机抽取一个数
int random(double* prob, const int n) {
  int i;
  double p = rand() / 1.0 / (RAND_MAX + 1);
  for (i = 0; i < n - 1; i++) {
    if (p <= prob[i])
      break;
    p -= prob[i];
  }
  return i;
}

// 根据转移矩阵和发射矩阵生成一串随机状态和符号
void randSeq(State* st, Result* res, const int n) {
  int i, ls, lr;
  srand((unsigned int) time(NULL));
  ls = random(init, nstate);
  lr = random(emission[ls], nresult);
  st[0] = state[ls];
  res[0] = result[lr];
  for (i = 1; i < n; i++) {
    ls = random(trans[ls], nstate);
    lr = random(emission[ls], nresult);
    st[i] = state[ls];
    res[i] = result[lr];
  }
}

void printState(State* st, const int n) {
  int i;
  for (i = 0; i < n; i++)
    printf("%c", st[i]);
  printf("\n");
}

void printResult(Result* res, const int n) {
  int i;
  for (i = 0; i < n; i++)
    printf("%c", res[i]);
  printf("\n");
}