原創: 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算法和后验解码的比较的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。