博客
关于我
强烈建议你试试无所不能的chatGPT,快点击我
Gibbs sampling
阅读量:5274 次
发布时间:2019-06-14

本文共 2802 字,大约阅读时间需要 9 分钟。

关于Gibbs sampling, 首先看一下上的解释:Gibbs sampling or Gibbs sampler is an  to generate a sequence of samples from the joint probability distribution of two or more . The purpose of such a sequence is to approximate the joint distribution, or to compute an  (such as an ).

说到Gibbs Sampling 就不得不说markov chain了。

Markov chain 是一组事件的集合,在这个集合中,事件是一个接一个发生的,并且下一个事件的发生,只由当前发生的事件决定。用数学符号表示就是:

A={ a1,a2 … ai, ai+1,… at }

P(ai+1| a1,a2,…ai) = P(ai+1| ai)

这里的ai不一定是一个数字,它有可能是一个向量,或者一个矩阵,例如我们比较感兴趣的问题里ai=(g, u, b)这里g表示基因的效应,u表示环境效应,b表示固定效应,假设我们研究的一个群体,g,u,b的联合分布用π(a)表示。事实上,我们研究QTL,就是要找到π(a),但是有时候π(a)并不是那么好找的,特别是我们要估计的a的参数的个数多于研究的个体数的时候。用一般的least square往往效果不是那么好。

解决方案:

用一种叫Markov chain Monte Carlo (MCMC)的方法产生Markov chain,产生的Markov chain{a1,a2 … ai, ai+1,… at }具有如下性质:当t 很大时,比如10000,那么at ~ π(a),这样的话如果我们产生一个markov chain:{a1,a2 … ai, ai+1,… a10000},那么我们取后面9000个样本的平均

  a_hat = (g_hat,u_hat,b_hat) = ∑a/ 9000 (i=1001,1002, … 10000)

这里g_hat, u_hat, b_hat 就是基因效应,环境效应,以及固定效应的估计值

MCMC有很多算法,其中比较流行的是Metropolis-Hastings Algorithm,Gibbs Sampling是Metropolis-Hastings Algorithm的一种特殊情况。MCMC算法的关键是两个函数:

1)    q(ai, ai+1),这个函数决定怎么基于ai得到ai+1

2)    α(ai, ai+1),这个函数决定得到的ai+1是否保留

目的是使得at的分布收敛于π(a)

Gibbs Sampling的算法:

一般来说我们通常不知道π(a),但我们可以得到p(g | u , b),p(u | g , b), p ( b | g, u )即三个变量的posterior distribution

Step1: 给g, u, b 赋初始值:(g0,u0,b0

Step2: 利用p (g | u0, b0) 产生g1

Step3: 利用p (u | g1, b0) 产生u1

Step4: 利用p (b | g1, u1) 产生b1

Step5: 重复step2~step5 这样我们就可以得到一个markov chain {a1,a2 … ai, ai+1,… at}

这里的q(ai, ai+1)= p(g | u , b)* p(u | g , b)* p ( b | g, u )

这里通俗点的解释一下。首先,什么是sampling。sampling就是以一定的概率分布,看发生什么事件。举一个例子。甲只能E:吃饭、学习、打球,时间T:上午、下午、晚上,天气W:晴朗、刮风、下雨。现在要一个sample,这个sample可以是:打球+下午+晴朗。。。 问题是我们不知道p(E,T,W),或者说,不知道三件事的联合分布。当然,如果知道的话,就没有必要用gibbs sampling了。但是,我们知道三件事的conditional distribution。也就是说,p(E|T,W),p(T|E,W),p(W|E,T)。现在要做的就是通过这三个已知的条件分布,再用gibbs sampling的方法,得到joint distribution。 具体方法。首先随便初始化一个组合,i.e. 学习+晚上+刮风,然后依条件概率改变其中的一个变量。具体说,假设我们知道晚上+刮风,我们给E生成一个变量,比如,学习-》吃饭。我们再依条件概率改下一个变量,根据学习+刮风,把晚上变成上午。类似地,把刮风变成刮风(当然可以变成相同的变量)。这样学习+晚上+刮风-》吃饭+上午+刮风。 同样的方法,得到一个序列,每个单元包含三个变量,也就是一个马尔可夫链。然后跳过初始的一定数量的单元(比如100个),然后隔一定的数量取一个单元(比如隔20个取1个)。这样sample到的单元,是逼近联合分布的。

=============================================

互动百科这么说的:

概率推理的通用方法,是Metropolis-Hastings算法的一个特例,因此也是Markov chain Monte Carlo算法的一种。

虽然它的通用性比较好,但导致了计算代价较高,所以在许多应用里,包括具有不完备信息的应用,都采用其它更为高效的方法。然而,理解这一方法有助于增进对统计推理问题的理解。

中心思想

由一个具有2个或更多变量的联合概率分布P(x1,x2,…,xn),生成一个样本序列{y1,y2,…,ym},用于逼近这一个联合分布,或计算一个积分(例如期望)。

适用于处理不完备信息,当联合分布不明确,而各个变量的条件分布已知的情况。

根据其他变量的当前值,依次对分布的每个变量生成一个实例。

随机过程

对一个随机过程,例如马尔可夫链过程,一般包括一个有限的状态集合 和 一个概率转移矩阵。假设这个过程各个各个状态都是可遍历的(ergodic),即转移矩阵中的元素值都大于0。为此,我们可以选择任意状态为初始态 Q0,计算转化N次后可能到达的状态 Qn 的概率。当N取值足够大时,可以计算得到这一过程最有可能的终态。

假设有一个变量集合X={X1,X2,……,Xn},P(X)为集合X的联合分布,0<p(x)<1。< p="">

我们将这些变量看做一个马尔科夫过程中的状态集,这一过程定义为:S=∏i=1~n

转载于:https://www.cnblogs.com/ywl925/archive/2013/03/19/2969134.html

你可能感兴趣的文章
CF219D Choosing Capital for Treeland
查看>>
杂七杂八的小笔记本
查看>>
51Nod1353 树
查看>>
CF1215E Marbles
查看>>
BZOJ2339 HNOI2011卡农(动态规划+组合数学)
查看>>
octave基本操作
查看>>
axure学习点
查看>>
WPF文本框只允许输入数字[转]
查看>>
dom4j 通用解析器,解析成List<Map<String,Object>>
查看>>
第一个项目--用bootstrap实现美工设计的首页
查看>>
使用XML传递数据
查看>>
TYVJ.1864.[Poetize I]守卫者的挑战(概率DP)
查看>>
0925 韩顺平java视频
查看>>
iOS-程序启动原理和UIApplication
查看>>
mysql 8.0 zip包安装
查看>>
awk 统计
查看>>
模板设计模式的应用
查看>>
实训第五天
查看>>
平台维护流程
查看>>
2012暑期川西旅游之总结
查看>>