什么是狄利克雷先验

机器算法验证 可能性 生物统计学 狄利克雷分布
2022-03-31 14:05:42

我正在做一些生物信息学研究,但我的背景是应用数学,我通常在其中的统计部分遇到问题。

基本上,我使用 bioStrings 包(bioConductor R 生物信息学套件)中的 R 函数创建了一个位置权重矩阵。

现在,读取函数签名看起来像

PWM(x, type = c("log2probratio", "prob"), prior.params = c(A=0.25, C=0.25, G=0.25, T=0.25))

对于type参数,详细信息部分说明了这一点

PWM 函数使用具有 Dirichlet 共轭的多项式模型来计算基 b 在位置 i 的估计概率。如参数部分所述,prior.params 提供 Dirichlet 先验中 DNAbase A、C、G 和 T 的参数。这些值导致碱基概率的位置独立初始估计为priorProbs = prior.params/sum(prior.params)和每个位置的碱基概率的后验(数据注入)估计为postProbs = (consensusMatrix(x) + prior.params)/(length(x) + sum(prior.params))type = "log2probratio"PWM=unitScale(log2(postProbs/priorProbs))type = "prob"PWM = unitScale(postProbs)

有人可以帮我理解这里的两种不同类型吗?如果位置频率矩阵的条目 有很多零(或小数据集),通常 PWM 将添加一个伪计数。PFM 的列将乘以 Dirichlet 分布。我想我真的只需要关于 Dirichlet 分布是什么以及它与多项模型的关系的帮助。

2个回答

让我试着回答你关于理解狄利克雷分布的最后一个问题,它与多项式的关系,我怀疑你真正想知道的是如何在应用环境中解释这一点,比如你的基因组学问题。

现在我将使用我对单倍群 SNP 的模糊回忆来解释所有这些,这可能与您的数据有些相似:

  • 所以假设我有这个数据集,灵感来自智人的这篇随机 NIH 论文,我需要识别与我认为存在但我不知道有多少的单倍群 N 的这个新颖的亚亚亚支相关的所有 SNP的 SNP (或者我猜,另一个群体遗传学术语是发现连锁不平衡)在这 4 个亚亚进化枝中的每一个中都有,以及这些 SNP 在这个 Y-DNA 样本中是什么。

15121 caacagcctt cataggctat gtcctcccgt gaggccaaat atcattctga ggggccacag

第15181章

第15241章

15301 tgcccttcat tattgcagcc ctagcagcac tccacctcct attcttgcac gaaacgggat

……等等……等等……`

  • 鉴于我们不知道连接在一起的 SNP 的数量会落入这个亚亚亚进化枝,或者在我即将发现的这个新的模糊基因组区域中出现 SNP 的概率,我将处理那些未知参数在贝叶斯范式下随机:

  • 我实际上会估计**与该亚亚亚进化相关联的 SNP 的数量**,因为我知道如果超过 1% 的群体在 DNA 序列的特定位置不携带相同的核苷酸, 那么这个变异可以归类为一个 SNP。

  • “关联的 SNP”是指一些未知数量的 SNP 组,例如,我们正在考虑与单倍群 N 的这个亚亚亚进化基因组区域相关的一个可能组是多巴胺能受体的 SNP 组x、y 和 z,另一组 - 用于血清素 5-HT2A 受体,即 SNP rs6311 和 rs6313 等)

  • 个采样核苷酸上观察到结果 SNP的预期次数(由参数表示,其中kk2iN

    X1,,Xk,xi(0,1)
    一个随机类别计数向量 ,
    i=1Nxi=1

由伪计数参数参数化α=(α1,,αk)

现在,一分钟的数学:

  1. 众所周知的概率分布是相关的,因为它在公共分布之间的关系图上清楚地说明,从我称之为“统计圣经”的Leemis(1986) 中采用,每个统计学家在研究生数学统计考试前都会看到生动的梦想学校,又名G.Casella 和 R.Berger 的“统计推理”第二版,2001 年。Cengage Learning。. 尽管即使在美国统计学家的地图的扩展版本中,也没有描绘狄利克雷和多项式,但这是我对它们在经典地图 #1 中的放置位置的大胆看法:

狄利克雷在哪里

此外,沿着这些箭头的某处将是另一个相关分布,即分类分布。

  1. 简而言之,如果它们在地图上非常接近并且来自同一个分布族,您可以导出或轻松将它们转换为另一个。其中一些是其他分布的推广,因此包括例如狄利克雷,它是对贝塔分布的推广,即狄利克雷将贝塔推广到多个维度

  2. 由于这个原因以及许多其他原因,Dirichlet 分布是多项分布的先验

  3. 现在回到我们的 SNP 问题:

我们将用于我们的问题的设置基于Pólya urn 模型 my-sub-sub-subclade-linked SNPs中显示的 4 字母核苷酸碱基的替换个字符串进行采样,其中对观察到的核苷酸碱基的每个 SNP 进行采样Nk p1,,pk

由于我们不知道有多少和哪些 SNP 属于“此子亚子进化枝的关联 SNP”,我们假设该子亚进化枝可能存在或不存在的那些未知 SNP 由参数,他和那些实际上可能是您引用的伪计数。这是一个很好的回应,并参考了这种 Dirichlet-Miltonomial 方法的一些问题。我强烈建议检查 Bioconductor 或该软件包文档,因为这些伪计数很容易成为转换不同矩阵同时集成非常不同分布的简单方法。αα1,,αk

  1. 现在我们在 Dirichlet-multinomial 模型中估计参数并更新它们:最终获得后验分布并估计属于我的新型子进化枝的链接 SNP 类型的数量用他们的概率来推断我们在子进化枝中看到这些群体的可能性有多大。α1+n1,,αk+nk

PS 我怀疑 Dirichlet 和 Multinomial 如此适用于基因组学并用于 R 中的一些生物信息学包的原因可能是由于这些类型的数据集以及传统模型(例如基于Hardy-Weinberg的模型)中的大量离散化原则,因为它们在数学上是某些二项式、Beta、多项式等设置类型的完美候选者,因为您实际上是在估计某些离散类别中某些计数的频率(尽管 Dirichlet 不要求参数是整数。α

PPS 从上面的一个非常简化的设置中,我省略了其他重要的事情,可以在本教程中找到,例如 Dirichlet 过程,特别是一些概率空间的划分步骤找到Θ

θk|Prior distribution H over component parameters, θkH

这个演示文稿中有一个很好的解释。

https://www.slideshare.net/g33ktalk/machine-learning-meetup-12182013

如果您愿意,您可以观看整个演示文稿(这是对 Dirichlet 分布的一个很好的解释),但我认为幻灯片会很快地传达这个概念。

幻灯片 32-35 解释了狄利克雷先验的数学过程。

幻灯片 50-60 显示了分发更新时发生的情况并显示了先前的情况。(视觉上比解释更容易)这可以理解总体思路

幻灯片 94-102 显示了更新发生时整个系统会发生什么。这与幻灯片 50-60 的概念相同,但会跟踪每次迭代发生的情况。