原创：hxj7
本文比较了viterbi算法求解最可能路径以及后验解码这两种不同的解码方法。
前文《序列比对（十）viterbi算法求解最可能路径》介绍了用viterbi算法求解最可能路径：在符号序列已知而状态序列未知时，最可能路径是：
image
本文将这两种方法比较了以下，看它们各自求解的路径差异是否显著。分两种情况：
一、如前面几篇文章一样，从公平骰子转为作弊骰子的概率是0.05。
效果如下：（其中Rolls一行是符号序列，也就是骰子投出的结果；Die一行是真实的骰子状态；Viterbi一行是viterbi算法求解出的最可能路径；PostDec一行是后验解码得出的路径）
image
二、将公平骰子转为作弊骰子的概率改为0.01。并将投骰子的次数增加到1000次。《生物序列分析》一书中说，此种情况下，viterbi求解的路径没有出现过'L'（即作弊骰子）。但是，笔者实验的结果是两种方法都可能出现'L'。效果如下：
image
具体代码如下：（以概率0.01，投骰子次数1000的情形为例写的代码）
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#define MIN_LOG_VALUE -99
//#define SAFE_EXP(x) ((x) < MIN_LOG_VALUE ? 0 : exp(x))
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.99, 0.01,
0.1, 0.9
};
const int nstate = 2;
const int nresult = 6;
double** fscore; // 前向算法的得分矩阵
double** bscore; // 后向算法的得分矩阵
double* scale; // 缩放因子向量
double logScaleSum;
State* rst; // 一串随机状态序列
Result* rres; // 一串随机符号序列
State* vst; // viterbi算法猜出来的状态序列
State* pst; // 后验解码得到的状态序列
struct Unit {
double v;
int *p;
int size;
};
typedef struct Unit* pUnit;
int random(double* prob, const int n);
void randSeq(State* st, Result* res, const int n);
int getResultIndex(Result r);
void printState(State* st, const int n);
void printResult(Result* res, const int n);
double forward(Result* res, const int n);
double backward(Result* res, const int n);
double** getPostProb(const int n);
void postDecode(double** prob, const int n);
void traceback(pUnit** a, const int l, const int i, State* st, const int m, int n);
void viterbi(Result* res, State* gst, const int n);
int main(void) {
int i;
int n = 1000;
double** postProb;
if ((rst = (State*) malloc(sizeof(State) * n)) == NULL || \
(rres = (Result*) malloc(sizeof(Result) * n)) == NULL || \
(scale = (double*) malloc(sizeof(double) * n)) == NULL || \
(fscore = (double**) malloc(sizeof(double*) * nstate)) == NULL || \
(bscore = (double**) malloc(sizeof(double*) * nstate)) == NULL || \
(vst = (Result*) malloc(sizeof(Result) * n)) == NULL) {
fputs("Error: out of space!\n", stderr);
exit(1);
}
for (i = 0; i < nstate; i++) {
if ((fscore[i] = (double*) malloc(sizeof(double) * n)) == NULL || \
(bscore[i] = (double*) malloc(sizeof(double) * n)) == NULL) {
fputs("Error: out of space!\n", stderr);
exit(1);
}
}
randSeq(rst, rres, n);
//printState(rst, n);
//printResult(rres, n);
forward(rres, n);
backward(rres, n);
postProb = getPostProb(n);
postDecode(postProb, n);
viterbi(rres, vst, n);
free(rst);
free(rres);
free(scale);
free(fscore);
free(bscore);
free(vst);
free(pst);
for (i = 0; i < nstate; i++)
free(postProb[i]);
free(postProb);
}
// 根据一个概率向量从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];
}
}
int getResultIndex(Result r) {
return r - result[0];
}
// 前向算法计算P(x)
double forward(Result* res, const int n) {
int i, l, k, idx;
double logpx;
// 缩放因子向量初始化
for (i = 0; i < n; i++)
scale[i] = 0;
// 计算第0列分值
idx = getResultIndex(res[0]);
for (l = 0; l < nstate; l++) {
fscore[l][0] = emission[l][idx] * init[l];
scale[0] += fscore[l][0];
}
for (l = 0; l < nstate; l++)
fscore[l][0] /= scale[0];
// 计算从第1列开始的各列分值
for (i = 1; i < n; i++) {
idx = getResultIndex(res[i]);
for (l = 0; l < nstate; l++) {
fscore[l][i] = 0;
for (k = 0; k < nstate; k++) {
fscore[l][i] += fscore[k][i - 1] * trans[k][l];
}
fscore[l][i] *= emission[l][idx];
scale[i] += fscore[l][i];
}
for (l = 0; l < nstate; l++)
fscore[l][i] /= scale[i];
}
// P(x) = product(scale)
// P(x)就是缩放因子向量所有元素的乘积
logpx = 0;
for (i = 0; i < n; i++)
logpx += log(scale[i]);
//printf("forward: logP(x) = %f\n", logpx);
logScaleSum = logpx;
/*
for (l = 0; l < nstate; l++) {
for (i = 0; i < n; i++)
printf("%f ", fscore[l][i]);
printf("\n");
}
*/
return exp(logpx);
}
// 后向算法计算P(x)
// backward算法中使用的缩放因子和forward中的一样
double backward(Result* res, const int n) {
int i, l, k, idx;
double tx, logpx;
// 计算最后一列分值
for (l = 0; l < nstate; l++)
bscore[l][n - 1] = 1 / scale[n - 1];
// 计算从第n - 2列开始的各列分值
for (i = n - 2; i >= 0; i--) {
idx = getResultIndex(res[i + 1]);
for (k = 0; k < nstate; k++) {
bscore[k][i] = 0;
for (l = 0; l < nstate; l++) {
bscore[k][i] += bscore[l][i + 1] * trans[k][l] * emission[l][idx];
}
}
for (l = 0; l < nstate; l++)
bscore[l][i] /= scale[i];
}
// 计算P(x)
tx = 0;
idx = getResultIndex(res[0]);
for (l = 0; l < nstate; l++)
tx += init[l] * emission[l][idx] * bscore[l][0];
logpx = log(tx) + logScaleSum;
//printf("backward: logP(x) = %f\n", logpx);
/*
for (l = 0; l < nstate; l++) {
for (i = 0; i < n; i++)
printf("%f ", bscore[l][i]);
printf("\n");
}
*/
return exp(logpx);
}
// 计算后验概率
double** getPostProb(const int n) {
int i, k;
double** postProb;
//double logdiff;
if ((postProb = (double**) malloc(sizeof(double*) * nstate)) == NULL) {
fputs("Error: out of space!\n", stderr);
exit(1);
}
for (k = 0; k < nstate; k++) {
if ((postProb[k] = (double*) malloc(sizeof(double) * n)) == NULL) {
fputs("Error: out of space!\n", stderr);
exit(1);
}
}
// 计算后验概率
for (i = 0; i < n; i++) {
for (k = 0; k < nstate; k++) {
postProb[k][i] = scale[i] * fscore[k][i] * bscore[k][i];
}
}
/*
printf("\n");
printf("Posterior Probabilities:\n");
for (k = 0; k < nstate; k++) {
for (i = 0; i < n; i++)
printf("%f ", postProb[k][i]);
printf("\n");
}
*/
return postProb;
}
void postDecode(double** prob, const int n) {
int i, k;
double maxCol;
int idx;
State* st;
if ((st = (State*) malloc(sizeof(State) * n)) == NULL) {
fputs("Error: out of space!\n", stderr);
exit(1);
}
for (i = 0; i < n; i++) {
idx = 0;
maxCol = prob[0][i];
for (k = 1; k < nstate; k++)
if (prob[k][i] > maxCol) {
maxCol = prob[k][i];
idx = k;
}
st[i] = state[idx];
}
/*
printf("\n");
printf("Posterior Decode:\n");
printState(st, n);
*/
pst = st;
}
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");
}
void traceback(pUnit** a, const int l, const int i, State* st, const int m, int n) {
int j, k;
int ll = 125; // 每行打印几个元素
int nl, nd;
pUnit pu = a[l][i];
if (! i) {
st[n] = state[l];
nl = m / ll;
nd = m % ll;
for (k = 0; k < nl; k++) {
printf("Rolls\t");
printResult(rres + k * ll, ll);
printf("Die\t");
printState(rst + k * ll, ll);
printf("Viterbi\t");
printState(st + k * ll, ll);
printf("PostDec\t");
printState(pst + k * ll, ll);
printf("\n");
}
if (nd > 0) {
printf("Rolls\t");
printResult(rres + k * ll, nd);
printf("Die\t");
printState(rst + k * ll, nd);
printf("Viterbi\t");
printState(st + k * ll, nd);
printf("PostDec\t");
printState(pst + k * ll, nd);
printf("\n");
}
printf("\n\n");
return;
}
st[n] = state[l];
for (j = 0, k = 0; j < nstate && k < pu->size; j++) {
if (pu->p[j]) {
traceback(a, j, i - 1, st, m, n - 1);
k++;
}
}
}
void viterbi(Result* res, State* gst, const int n) {
double maxCol;
double* tm;
int i, j, k, l;
int idx;
pUnit** aUnit; // 得分矩阵
double* loginit; // 每个元素都取log后的初始向量
double** logem; // 每个元素都取log后的发射矩阵
double** logtrans; // 每个元素都取log后的转移矩阵
double v0 = 0; // v0(0)的log值
// 初始化
if ((aUnit = (pUnit**) malloc(sizeof(pUnit*) * nstate)) == NULL || \
(loginit = (double*) malloc(sizeof(double) * nstate)) == NULL || \
(logem = (double**) malloc(sizeof(double*) * nstate)) == NULL || \
(logtrans = (double**) malloc(sizeof(double*) * nstate)) == NULL) {
fputs("Error: out of space!\n", stderr);
exit(1);
}
for (i = 0; i < nstate; i++) {
if ((aUnit[i] = (pUnit*) malloc(sizeof(pUnit) * n)) == NULL || \
(logem[i] = (double*) malloc(sizeof(double) * nresult)) == NULL || \
(logtrans[i] = (double*) malloc(sizeof(double) * nstate)) == NULL) {
fputs("Error: out of space!\n", stderr);
exit(1);
}
for (j = 0; j < n; j++) {
if ((aUnit[i][j] = (pUnit) malloc(sizeof(struct Unit))) == NULL || \
(aUnit[i][j]->p = (int*) malloc(sizeof(int) * nstate)) == NULL) {
fputs("Error: out of space!\n", stderr);
exit(1);
}
for (k = 0; k < nstate; k++)
aUnit[i][j]->p[k] = 0;
aUnit[i][j]->size = 0;
}
}
if ((tm = (double*) malloc(sizeof(double) * nstate)) == NULL) {
fputs("Error: out of space!\n", stderr);
exit(1);
}
// 初始向量取log值
for (i = 0; i < nstate; i++)
loginit[i] = init[i] == 0 ? MIN_LOG_VALUE : log(init[i]);
// 发射矩阵取log值
for (i = 0; i < nstate; i++)
for (j = 0; j < nresult; j++)
logem[i][j] = emission[i][j] == 0 ? MIN_LOG_VALUE : log(emission[i][j]);
// 转移矩阵取log值
for (i = 0; i < nstate; i++)
for (j = 0; j < nstate; j++)
logtrans[i][j] = trans[i][j] == 0 ? MIN_LOG_VALUE : log(trans[i][j]);
// 动态规划计算得分矩阵
// 首先计算第0列，因为第0列的值和vk(0)有关
// v0(0) = 1, vk(0) = 0 for k>0
idx = getResultIndex(res[0]);
for (l = 0; l < nstate; l++)
aUnit[l][0]->v = v0 + loginit[l] + logem[l][idx];
// 计算从第1列开始的各列
for (i = 1; i < n; i++) {
idx = getResultIndex(res[i]);
for (l = 0; l < nstate; l++) {
maxCol = tm[0] = aUnit[0][i - 1]->v + logtrans[0][l];
for (k = 1; k < nstate; k++) {
tm[k] = aUnit[k][i - 1]->v + logtrans[k][l];
if (tm[k] > maxCol)
maxCol = tm[k];
}
aUnit[l][i]->v = maxCol + logem[l][idx];
for (k = 0; k < nstate; k++)
if (tm[k] == maxCol) {
aUnit[l][i]->p[k] = 1;
aUnit[l][i]->size++;
}
}
}
/*
// 打印得分矩阵
for (l = 0; l < nstate; l++) {
for (i = 0; i < n; i++)
printf("%f ", aUnit[l][i]->v);
printf("\n");
}
*/
maxCol = aUnit[0][n - 1]->v;
for (l = 1; l < nstate; l++)
if (aUnit[l][n - 1]->v > maxCol)
maxCol = aUnit[l][n - 1]->v;
for (l = 0; l < nstate; l++)
if (aUnit[l][n - 1]->v == maxCol) {
traceback(aUnit, l, n - 1, gst, n, n - 1);
}
}
（公众号：生信了）
网友评论