原創：hxj7
本文介紹如何計算狀態的後驗概率。
的概率。很明顯，此概率為一後驗概率。
要計算上述後驗概率，可以經過以下推導：
其中：
根據公式(1),(4),(5),(6)，可以重新計算後驗概率：
據公式(7)，後驗概率計算就簡單多了。可以利用前文**，稍加增改即可。執行效果如下：
state state = ; // 所有的可能狀態
result result = ; // 所有的可能符號
double init = ; // 初始狀態的概率向量
double emission[6] = ;
double trans[2] = ;
const int nstate = 2;
const int nresult = 6;
double** fscore; // 前向演算法的得分矩陣
double** bscore; // 後向演算法的得分矩陣
double* scale; // 縮放因子向量
double logscalesum;
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);
int main(void)
for (i = 0; i < nstate; i++)
} randseq(rst, rres, n);
printstate(rst, n);
printresult(rres, n);
forward(rres, n);
backward(rres, n);
postprob = getpostprob(n);
free(rst);
free(rres);
free(scale);
free(fscore);
free(bscore);
for (i = 0; i < nstate; i++)
free(postprob[i]);
free(postprob);
}// 根據乙個概率向量從0到n-1隨機抽取乙個數
int random(double* prob, const int n)
return i;
}// 根據轉移矩陣和發射矩陣生成一串隨機狀態和符號
void randseq(state* st, result* res, const int n)
}int getresultindex(result r)
// 前向演算法計算p(x)
double forward(result* res, const int n)
for (l = 0; l < nstate; l++)
fscore[l][0] /= scale[0];
// 計算從第1列開始的各列分值
for (i = 1; i < n; i++)
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++)
*/ return exp(logpx);
}// 後向演算法計算p(x)
// backward演算法中使用的縮放因子和forward中的一樣
double backward(result* res, const int n)
}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++)
*/ return exp(logpx);
}// 計算後驗概率
double** getpostprob(const int n)
for (k = 0; k < nstate; k++)
} // 計算後驗概率
for (i = 0; i < n; i++)
} printf("n");
printf("posterior probabilities:n");
for (k = 0; k < nstate; k++)
return postprob;
}void printstate(state* st, const int n)
void printresult(result* res, const int n)
序列比對 levenshtein distance
一 兩個字串a和b，將a變為b需要經過最少幾次編輯 包括插入 刪除 替換 解決問題，俄羅斯人levenshtein設計了乙個演算法，名為levenshtein距離，具體如下 levenshtein的具體演算法是 1.獲取字串 的長度為n，獲取字串 的長度為m，如果m 0，返回n 如果n 0，返回m。...
生物序列比對與計算智慧型
生物資訊學的首要任務之一就是從資料庫中搜尋同源序列，尋找保守的序列模式。而序列比對是最常用的方法，對於發現生物序列中的功能 結構和進化的資訊具有非常重要的意義。下面給出序列比對的定義 定義 序列比對問題可以表示為乙個五元組 msa s a f 其中 1 為序列比對的符號集 表示空位 gap 表示基本...
Michael Schatz 序列比對課程
michael schatz cold spring harbor laboratory 最近在研究 bwa mem 序列比對演算法，直接去看 看不懂，就3頁，太精簡了，好多背景知識都不了解。通過google，發現了一系列序列比對的課程，講得真是太好了，真是撿到寶了！如果你想搞透生物資訊的基石 序列...