viterbi算法_序列比对(十四)——viterbi算法和后验解码的比较

原创: hxj7

本文比较了viterbi算法求解最可能路径以及后验解码这两种不同的解码方法。

前文《序列比对（十）viterbi算法求解最可能路径》介绍了用viterbi算法求解最可能路径：在符号序列已知而状态序列未知时，最可能路径是：

本文将这两种方法比较了一下，看它们各自求解的路径差异是否显著。分两种情况：

一、如前面几篇文章一样，从公平骰子转为作弊骰子的概率是0.05。
效果如下：（其中Rolls一行是符号序列，也就是骰子投出的结果；Die一行是真实的骰子状态；Viterbi一行是viterbi算法求解出的最可能路径；PostDec一行是后验解码得出的路径）

二、将公平骰子转为作弊骰子的概率改为0.01。并将投骰子的次数增加到1000次。《生物序列分析》一书中说，此种情况下，viterbi求解的路径没有出现过'L'（即作弊骰子）。但是，笔者实验的结果是两种方法都可能出现'L'。效果如下：

具体代码如下：（以概率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) = %fn", 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) = %fn", 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("Rollst");printResult(rres + k * ll, ll);printf("Diet");printState(rst + k * ll, ll);printf("Viterbit");printState(st + k * ll, ll);printf("PostDect");printState(pst + k * ll, ll);printf("n");}if (nd > 0) {printf("Rollst");printResult(rres + k * ll, nd);printf("Diet");printState(rst + k * ll, nd);printf("Viterbit");printState(st + k * ll, nd);printf("PostDect");printState(pst + k * ll, nd);printf("n"); }printf("nn");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>0idx = 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);}
}

（公众号：生信了）

viterbi算法_序列比对(十四)——viterbi算法和后验解码的比较相关推荐

中希尔排序例题代码_超全面分析十大排序算法
点击上方"零一视界",选择"星标"公众号资源干货,第一时间送达作者 | 不该相遇在秋天责编 | 程序员小吴前言本文全长 14237 字,配有 70 张 ...
学术会议墙报_中国化学会第十四届全国电分析化学学术会议在南京顺利召开
2020年11月26-29日,由中国化学会分析化学学科委员会和中国仪器仪表学会分析仪器分会主办,东南大学化学化工学院和南京大学生命分析化学国家重点实验室承办的中国化学会第十四届全国电分析化学学术会议在 ...
（十四）算法设计思想之“贪心算法”
算法设计思想之"贪心算法" 贪心算法是什么 LeetCode:455.分饼干 LeetCode:122.买卖股票的最佳时机II 思考题贪心算法是什么贪心算法是算法设计中的一种方 ...
数据结构与算法之美（十四）算法思想——贪心算法
目录贪心算法介绍贪心算法例子 1. 背包 2. 分糖果 3. 钱币找零 4. 区间覆盖 5. 区间覆盖的延伸:任务调度.教师排课贪心算法经典应用 1. 霍夫曼编码 2. 最小生成树算法 3. 最 ...
病虫害模型算法_基于深度学习的目标检测算法综述
sigai 基于深度学习的目标检测算法综述导言目标检测的任务是找出图像中所有感兴趣的目标(物体),确定它们的位置和大小,是机器视觉领域的核心问题之一.由于各类物体有不同的外观,形状,姿态,加上成像 ...
hash算法_到底什么是Hash？Hash算法的原理和实际应用讲解
提到hash,相信大多数同学都不会陌生,之前很火现在也依旧很火的技术区块链背后的底层原理之一就是hash,下面就从hash算法的原理和实际应用等几个角度,对hash算法进行一个讲解. 1.什么是Has ...
关联规则挖掘算法_基于Apriori关联规则的协同过滤算法
Apriori 算法 apriori关联规则算法的原理设计较为简单,著名的"啤酒和尿布"说的就是Apriori算法,通俗来讲apriori旨在寻找频繁项集,以帮助商家将消费者有可能 ...
蝴蝶优化算法_腾讯机智团队分享--AllReduce算法的前世今生
从事分布式深度学习相关工作的同学,应该都频繁地用到了AllReduce(规约)操作. 图1 AllReduce的示意图但是对于训练框架中集成的AllReduce相关操作,其背后实现的原理是什么? 除 ...
数据与广告系列二十四：效果广告后定向时代如何逆流而上
作者·黄崇远『数据虫巢』全文共4338字题图ssyer.com " 在效果广告的发展历程中,当前已经处于后定向时代,或者说是弱定向时代,我们是应该顺应潮流还是应该逆流而上?" ...

viterbi算法_序列比对(十四)——viterbi算法和后验解码的比较

viterbi算法_序列比对(十四)——viterbi算法和后验解码的比较相关推荐

最新文章

热门文章