序列每天从0开始_序列比对(十一)——计算符号序列的全概率

3b304757986189285d23965e1b7ddcad.png
前文介绍了在知道符号序列后用viterbi算法求解最可能路径。本文介绍了如何使用前向算法和后向算法计算符号序列的全概率。

如果一个符号序列中每个符号所对应的状态是已知的,那么这个符号序列出现的概率是容易计算的:

b86039f365a6ed14ba2bdca5fc9ccef2.png

但是,如果一个符号序列中每个符号所对应的状态未知时,该怎么求取这条序列的概率呢?我们知道:

dd7491c1b40841161490159b7f0866be.png

如果我们用穷举法求出所有的P(x,π)是不现实的,因为随着序列长度的增长,所有可能的路径的数目是指数增长的。这个时候,我们可以再次借助动态规划来求取。

有两种方法,前向法后向法。二者的区别前向法是从序列头部开始计算,逐步向序列尾部推进;而后向法是从序列尾部开始计算,逐步向序列头部推进。

前向法

定义:

1cac482828c9a0f99cfbfdb7516083b4.png
图片引自《生物序列分析》

1433e576891e9fc1478ca9bb0154a0a9.png
图片引自《生物序列分析》

后向法

1473d7fd6fc9486990201442a0f3d650.png
图片引自《生物序列分析》


解决下溢的问题

与《序列比对(十)viterbi算法求解最可能路径》一文中的viterbi算法相似,前向法和后向法也都涉及到下溢的问题。由于递归公式中涉及到加法,所以不能像《序列比对(十)viterbi算法求解最可能路径》中简单使用log变换。《生物序列分析》一书中给出了两种解决方法:

一是近似的log变换

3d1a85bb32ab289624b8e53091306de7.png
图片引自《生物序列分析》

二是使用一组缩放因子

3f86769e3818c429fd5e7a7abc3f029b.png
图片引自《生物序列分析》

实现代码和效果

下面的代码首先随机生成一个状态序列和相应的符号序列,然后根据前向法和后向法来计算符号序列的全概率。本文采用缩放因子来解决下溢的潜在问题。这样做有一个好处,就是此时所有缩放因子的乘积就等于P(x)。

效果如下:

76b3be84c6f54b6be10faa6698c2cd3d.png

具体代码:

#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");  
}

(公众号:生信了)

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/news/304358.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

SQL 2005 使用row_number来分页

今天研究了一下row_number,用它来返回特定行的记录感觉是非常方便的&#xff0c;所以就做了个分页的存储过程&#xff0c;但不知道性能较之top和游标之类的那个好 代码 createprocedure[dbo].[proc_TestPage]--表名 tablenamenvarchar(255), --排序字段 sortcolumnnvarchar(255…

债务大爆发,中国30%家庭不堪一击!

债务大爆发&#xff0c;30%中国家庭“不堪一击”&#xff01;从2007年到2016年&#xff0c;中国家庭的债务率翻了一倍多。已经有超过1/3的家庭属于高负债家庭。前不久&#xff0c;深圳中兴网信科技有限公司的一研发组主管欧某&#xff0c;以最决绝的方式&#xff0c;从中兴通迅…

腾讯35k招.NET Core开发,深扒这些技术要求 真的很难吗?

3月草长莺飞&#xff0c;3月招聘满天飞&#xff0c;各种高薪招聘更是心里种草&#xff0c;前几天分享了腾讯牛年35k的.NET Core招聘需求&#xff0c;分享了一波资料深受好评&#xff0c;本着再接再厉的精神&#xff0c;本文继续为大家上干货&#xff0c;补齐腾讯的各种要求。新…

【直观理解】为什么梯度的负方向是局部下降最快的方向?

推荐阅读时间&#xff1a;8min~15min主要内容&#xff1a;为什么梯度的负方向是局部下降最快的方向&#xff1f;刚接触梯度下降这个概念的时候&#xff0c;是在学习机器学习算法的时候&#xff0c;很多训练算法用的就是梯度下降&#xff0c;然后资料和老师们也说朝着梯度的反方…

紫光物联linux登录账号,紫光展锐打造操作系统生态,赋能万物互联智能时代

本周&#xff0c;以“象由芯生科技服务人民”为主题的2020紫光展锐市场峰会重磅开启&#xff0c;广大生态合作伙伴共聚一堂&#xff0c;共话数字世界新未来。在今天举办的“操作系统OS研讨会”上&#xff0c;来自紫光展锐工程一线的架构师带来了一场整个操作系统领域的饕餮盛宴…

一行命令搭建内部的管道

在上一篇《边缘计算k8s集群之SuperEdge》文章中&#xff0c;笔者基于ECK搭建了边缘集群并添加了节点。通过边缘集群&#xff0c;我们可以很方便的管理各个地域的节点&#xff0c;本地、各云厂商的机房、客户所在地、海外的都可以。在本篇内容&#xff0c;我们将讲述如何使用ips…

ArchiMate - 发布【企业架构语言ArchiMate v0.5.pdf】

在《年度总结和计划&#xff1a;去年4个1&#xff0c;今年5个1》中说过今年我准备在项目组引入1个架构语言&#xff08;ArchiMate&#xff09;&#xff0c;为了便于大家学习&#xff0c;我把一些内容集成一本电子书&#xff0c;目前发布0.5版本&#xff0c;后续还会不断更新&am…

那些有趣/用的 Python 库

图片处理pip install pillowfrom PIL import Imageimport numpy as npa np.array(Image.open(test.jpg))b [255,255,255] - aim Image.fromarray(b.astype(uint8))im.save(new.jpg)youtube-dl下载国外视频pip install youtube-dl #直接安装youtube-dlpip install -U youtube…

linux系统刷分辨率,Linux下设置其分辨率及刷新率

行频&#xff1a;行频又称为“水平扫描频率”&#xff0c;指电子枪每秒在荧光屏上扫过的水平线的数量&#xff0c;其值等于“场频 垂直分辨率1.04”&#xff0c;单位为KHz(千赫兹)。行频是一个综合分辨率和场频的参数&#xff0c;该值越大&#xff0c;显示器可以提供的分辨率越…

.NET 5 部署在docker上运行

1、创建站点创建一个ASP.NET Core Web应用程序&#xff0c;选中启用Docker支持。自动帮我们创建一个Dockerfile文件。2、编写Dockerfile文件dockerfile是一个文件格式的配置文件&#xff0c;用户可以使用dockerfile来快速构建自定义的镜像。由一行行命令语句组成&#xff0c;并…

Nexus:一站式私有仓库管理(NuGet、Maven、npm、Docker)

我们在日常开发中经常需要使用到私有仓库&#xff0c;比如 dotNET 中的 NuGet、Java 中的 Maven、前端的 npm&#xff0c;还有 Docker 镜像&#xff0c;每一个私有仓库各自管理&#xff0c;维护起来比较麻烦&#xff0c;而 Nexus 可以将其统一起来。本文将介绍 Nexus 的安装以及…

众里寻 Bug 千百度,蓦然回首,它却在隔壁老张处……

程序员与 Bug 是一对矛盾的存在&#xff0c;程序员既要在解决 Bug 中获得成就感&#xff0c;同时也讨厌 Bug 本身的存在。“程序不息&#xff0c;Bug 不止”&#xff0c;程序员在与 Bug 的斗争中&#xff0c;也有很多有趣的事情发生&#xff0c;我们整理了一些程序员在调试 Bug…

Blazor WASM 实现人民币大写转换器

点击上方蓝字关注“汪宇杰博客”导语.NET 5 正式发布已经有一段时间了&#xff0c;其中 Blazor 技术是该版本的亮点之一。作为微软技术的被坑者&#xff0c;年少的我曾经以为 SilverLight 能血虐 Flash&#xff0c;Zune 能团灭 iPod&#xff0c;WP 能吊打 iPhone&#xff0c;UW…

金山安全实验室公布中国互联网六大类钓鱼网站

金山安全实验室公布中国互联网六大类钓鱼网站金山安全实验室研究人员对中国大陆钓鱼网站的普遍特征进行分析&#xff0c;发现以下六个领域最容易被钓鱼网站***&#xff1a;1.QQ十年庆典、QQ抽奖、腾讯活动&#xff1b;2.证券、股票分析、黑庄、理财专家等财经领域&#xff1b;3…

程序员过关斩将--领导说我的类的职责不单一

“为什么类的职责要单一化&#xff1f;“类的职责单一化很容易吗&#xff1f;首先&#xff0c;我要提醒一下看到这篇文章的同学&#xff0c;我认为保证类&#xff08;一定是类吗&#xff1f;&#xff09;的单一职责并不容易软件开发过程中&#xff0c;自古就流传着几大规则&…

从概念到案例,机器学习应该掌握的20个知识点

随着科技的发展&#xff0c;计算机对人类的生产活动和社会活动产生了极为重要的影响&#xff0c;同时以强大的生命力飞速发展着。目前计算机正广泛用于社会各个领域&#xff0c;并朝着微型化、网络化、智能化和巨型化的方向前进。说到智能化&#xff0c;大家最先想到的应该就是…

闲来没事写个记事本玩玩!!!

这两天工作压力大&#xff0c;还好今天见着太阳了&#xff0c;这会儿没事写了个记事本&#xff0c;功能单一&#xff0c;适合练手&#xff0c;可能对新手会有所帮助&#xff0c;贴上来看看吧&#xff0c; 说到实现 记事本&#xff0c;我们应该选择什么样的控件呢&#xff0c;Te…

LG将授权webOS给其他电视厂商使用

喜欢就关注我们吧&#xff01;LG 将向其他公司提供 webOS。根据 LG 发布的公告&#xff0c;其自家电视机搭载的专有系统 webOS 将会授权给其他的外部电视厂商使用。被授权使用 webOS 的电视厂商还会获得来自 LG 的 Magic Motion 遥控器&#xff0c;此外&#xff0c;系统的语音控…

数据之美,堪比好莱坞大片!

看完下面的几张图&#xff0c;你就知道自己有多无知了。堪称是好莱坞大片啊&#xff01;1城市3D空间通过2D瓦片图层的3D化&#xff0c;能够在经度维度、量级、时间多个维度上真实还原城市3D空间。例子中为模拟的轨迹数据和旧金山食物供应商分布。2GPS轨迹分布以三种不同的方式描…

HSRP的配置问题

HSRP的配置问题<?xml:namespace prefix o ns "urn:schemas-microsoft-com:office:office" />实验目的&#xff1a;理解和掌握路由热备份的配置步骤和原理实现网关的冗余功能实验环境&#xff1a;如下图所示<?xml:namespace prefix v ns "urn:sch…