序列比对(十二)——计算后验概率

最新推荐文章于 2024-04-10 21:23:25 发布

生信了（公众号同名）

最新推荐文章于 2024-04-10 21:23:25 发布

阅读量1.8k

点赞数

分类专栏： # 序列算法文章标签：二代测序生物信息序列比对算法

本文链接：https://blog.csdn.net/biocity/article/details/96830369

版权

本文详细阐述了如何在已知整条符号序列的情况下，通过数学推导和算法来计算序列中每个位置对应状态的后验概率。介绍了利用前向算法和后向算法的全概率，并给出公式及代码实现。

摘要由CSDN通过智能技术生成

原创: hxj7

本文介绍如何计算状态的后验概率。

前文《序列比对(十一)——计算符号序列的全概率》介绍了如何使用前向算法和后向算法计算符号序列的全概率。但是很多情况下我们也想了解在整条符号序列已知的情况下，某一位置符号所对应的状态的概率。也就是说要计算
在这里插入图片描述
的概率。很明显，此概率为一后验概率。
要计算上述后验概率，可以经过以下推导：

其中：

根据公式(1),(4),(5),(6)，可以重新计算后验概率：
在这里插入图片描述

据公式(7)，后验概率计算就简单多了。可以利用前文代码，稍加增改即可。运行效果如下：

**具体代码如下：

#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;

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) {
   
  int i;
  int n = 5;
  State* rst;   // 一串随机状态序列
  Result* rres;  // 一串随机符号序列
  double** postProb;
  if ((rst = (State