前文介紹了在知道符號序列后用viterbi算法求解最可能路徑。本文介紹了如何使用前向算法和后向算法計算符號序列的全概率。

如果一個符號序列中每個符號所對應的狀態是已知的，那么這個符號序列出現的概率是容易計算的：

但是，如果一個符號序列中每個符號所對應的狀態未知時，該怎么求取這條序列的概率呢？我們知道：

如果我們用窮舉法求出所有的P(x,π)是不現實的，因為隨著序列長度的增長，所有可能的路徑的數目是指數增長的。這個時候，我們可以再次借助動態規劃來求取。

有兩種方法，前向法和后向法。二者的區別是前向法是從序列頭部開始計算，逐步向序列尾部推進；而后向法是從序列尾部開始計算，逐步向序列頭部推進。

前向法

定義：

圖片引自《生物序列分析》

圖片引自《生物序列分析》

后向法

圖片引自《生物序列分析》

解決下溢的問題

與《序列比對（十）viterbi算法求解最可能路徑》一文中的viterbi算法相似，前向法和后向法也都涉及到下溢的問題。由于遞歸公式中涉及到加法，所以不能像《序列比對（十）viterbi算法求解最可能路徑》中簡單使用log變換。《生物序列分析》一書中給出了兩種解決方法：

一是近似的log變換

圖片引自《生物序列分析》

二是使用一組縮放因子

圖片引自《生物序列分析》

實現代碼和效果

下面的代碼首先隨機生成一個狀態序列和相應的符號序列，然后根據前向法和后向法來計算符號序列的全概率。本文采用縮放因子來解決下溢的潛在問題。這樣做有一個好處，就是此時所有縮放因子的乘積就等于P(x)。

效果如下：

具體代碼：

#include <stdio.h> #include <stdlib.h> #include <time.h> #include <math.h> //#define MIN_LOG_VALUE -15 //#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.95, 0.05,0.1, 0.9 }; const int nstate = 2; const int nresult = 6;State* rst; // 一串隨機狀態序列 Result* rres; // 一串隨機符號序列 double** fscore; // 前向算法的得分矩陣 double** bscore; // 后向算法的得分矩陣 double* scale; // 縮放因子向量 double logScaleSum;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 printResult(Result* res, const int n); double forward(const int n); double backward(const int n);int main(void) {int i;int n = 300;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) {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);printResult(rres, n);forward(n);backward(n);free(rst);free(rres);free(scale);free(fscore);free(bscore); }// 根據一個概率向量從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(const int n) {int i, l, k, idx;double logpx;// 縮放因子向量初始化for (i = 0; i < n; i++)scale[i] = 0;// 計算第0列分值idx = getResultIndex(rres[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(rres[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(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(rres[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(rres[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); }void printResult(Result* res, const int n) {int i;for (i = 0; i < n; i++)printf("%c", res[i]);printf("n"); }

（公眾號：生信了）

總結

以上是生活随笔為你收集整理的序列每天从0开始_序列比对(十一)——计算符号序列的全概率的全部內容，希望文章能夠幫你解決所遇到的問題。

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