随机模拟-Monte Carlo积分及采样(详述直接采样、接受-拒绝采样、重要性采样)

转载请注明出处 http://www.jianshu.com/p/3d30070932a8
作者:@贰拾贰画生


1. Monte Carlo 积分

蒙特卡洛方法的思想很简单,就是用随机投点法来模拟不规则图形的面积。
比如在1*1的矩形中,有一个不规则的图形,我们想要直接计算该图形的面积很困难,那怎么办呢?我们可以拿N个点,随机抛在1*1的矩形框中,数一下落入该不规则图形中的点的个数count,那么该不规则图形的面积就可以用count/N近似。

除了求面积,蒙特卡洛方法还有什么应用呢?

求积分。

有函数f(x),它在区间[0, a]的积分就是x=0, x=a以及f(x)和x轴围成区域的面积,我们将[0, a]划分n份,那么积分可以写成:



我们采用蒙特卡洛采样的思想,右边把极限去掉,只取n个点,那么剩下的部分岂不就是f(x)的均值E[f(x)]?没错,就是n个点的f(x)的均值。

所以,对于f(x)在[a, b]的积分,我们可以通过求f(x)的均值来模拟。

举个小例子,我们求这个积分:


我们随机取[0,1]上的N个点,然后求其均值就是模拟该积分。
MATLAB代码

N_samples = 100000;
f_N = zeros(N_samples, 1);
Z_N = rand(N_samples, 1);
 
for i = 1 : N_samples
    f_N (i, 1) = atan(Z_N(i, 1)) / (Z_N(i, 1) * Z_N(i, 1) + Z_N(i, 1) * sin(Z_N(i, 1)));
end
 
I_cap = mean(f_N)

Z_N保存随机取的N_samples = 100000个点的x值,f_N保存f(x),然后求f_N的均值,得到的就是积分I。


2. 直接采样

直接采样的思想是,通过对均匀分布采样,实现对任意分布的采样。因为均匀分布采样好猜,我们想要的分布采样不好采,那就采取一定的策略通过简单采取求复杂采样。
假设y服从某项分布p(y),其累积分布函数CDF为h(y),有样本z~Uniform(0,1),我们令 z = h(y),即 y = h(z)^(-1),结果y即为对分布p(y)的采样。

直接采样

举个例子:



直接采样有一个问题,就是上图下方的那两个问号。

3. 接收-拒绝采样

接受-拒绝采样的思想是,对于分布p(z),很难通过直接采样进行采样,但是我们有q(z),可以通过直接采样或其他采样进行采样,那么我们怎么可以通过q(z)采样得到p(z)的采样呢?
先看下图:

接收-拒绝采样

红色的是p(z), 蓝色的是q(z),我们对q(z)乘一个参数k,让k能正好包住p(z),那么对于每一个从q(z)得到的样本z0,我们有一定的概率接受它,概率的大小就是p(z0) / kq(z0)。很容易就能看出来,在p(z)和kq(z)相切的地方的采样,接受率就是1。那么有人问了,接受率能计算出来,但是我们对于一个样本z0,到底怎么判断是接受还是不接受啊?我们有u~Uniform[0,1],对于每一个样本z0,我们一个u0,如果u0 <= p(z0) / kq(z0),我们就接受,否则就拒绝。重复此过程,得到的样本就服从分布p(z)。

接受-拒绝采样算法

当然,q(z)的选取要有一定规则:q(z)与p(z)外形要相近,q(z)采样方便。

现在有一个问题,k怎么求?

其实也很简单,看图有 kq(z) >= p(z),那么k >= p(z) / q(z),我们求p(z) / q(z)的最大值,即为k。

举个例子,对截断正态分布的接受-拒绝采样。
截断正态分布的意思就是对于正态分布N(a, b) x属于[0, 4],其在[0, 4]上的积分为1,而不是在负无穷到正无穷的积分为1。截断正态分布不是正态分布,所以,我们知道截断正态分布的概率密度函数。

维基上对截断正态分布的定义:


截断正态分布

分子的小fai是标准正态分布的概率密度函数,分母上的大fai是标准正态分布的累积分布函数。

我们有p(z)服从N(1, 1), I(0 <= x <= 4),令q(z)~U[0, 4],根据上图的公式,p(z)就有了,q(z) = 1/4,所以k = max(p(z) / q(z)),在z = 1(均值)的时候p(z) / q(z)取最大,所以得到k:


这个例子比较巧,分母没有z,我们可以直接判断在z=1时,p(z)/q(z)取最大,如果p(z), q(z)都有z,那么要通过求导的方式求k了。

4. 重要性采样求均值

对于重要性采样,我之前是有个疑问的,既然是采样,为何最后没有得到样本,反而去求均值去了?其他很多介绍重要性采样的文章都没有讲明白这一点,其实重要性采样与接受-拒绝采样有异曲同工之妙。接受拒绝采样时通过接受拒绝的方式对通过q(z)得到的样本进行筛选使得最后得到的样本服从分布p(z),每个接受的样本没有高低贵贱之分,一视同仁。而重要性采样的思想是,对于通过q(z)得到的样本,我全部接受!全部接受的话会有一个问题,那就是最后样本点分布不能服从p(z)分布,为了矫正这个样本全部接受带来的偏差,我们给每个样本附一个重要性权重,比如,对于p(z0)/q(z0)=1的样本,给的权重大一点,p(z0)/q(z0)=0.1的样本,权重小一点。但是这个权重怎么算呢?哎,我们发现这个p(z0)/q(z0)不就是一个很好的权重标识吗?

重要性采样取得到的是带有重要性权重的服从q(z)分布的样本,这个权重乘以样本之后的结果其实就是服从p(z)分布的。对于这种比较特殊的样本,我们怎么用呢?我们一般用来求均值,所以,这一小节的标题是重要性采样求均值。如果只写重要性采样容易让人产生误解。

重要性采样

应当注意的是,图中p(z)与f(z)的关系,p(z)是一种分布,是相对于z轴的采样点而言的,比如在红色的两个驼峰处,z的取点比较多,在其他地方z的取点就比较少,这叫样本分布服从p(z)。对于f(z)是一种映射关系,将z值映射到其他维度。比如我们熟悉的y = f(x),将x映射到y。我们所说的求均值就是求f(z)的均值。这一点一定要想明白。

5. 马尔科夫蒙特卡洛采样 & Gibbs采样

这两部分,很多人已经写了不错的文章,我就不写了。
推荐一个博文MCMC(Markov Chain Monte Carlo) and Gibbs Sampling,介绍的已经很清楚了。

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 203,937评论 6 478
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 85,503评论 2 381
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 150,712评论 0 337
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,668评论 1 276
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,677评论 5 366
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,601评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,975评论 3 396
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,637评论 0 258
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,881评论 1 298
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,621评论 2 321
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,710评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,387评论 4 319
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,971评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,947评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,189评论 1 260
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 44,805评论 2 349
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,449评论 2 342

推荐阅读更多精彩内容

  • 查看原文 1 简介 Deep Learning最简单的一种方法是利用人工神经网络的特点,人工神经网络(ANN)本身...
    JinkeyAI阅读 6,732评论 0 4
  • Spring Cloud为开发人员提供了快速构建分布式系统中一些常见模式的工具(例如配置管理,服务发现,断路器,智...
    卡卡罗2017阅读 134,591评论 18 139
  • 今天利用早起1小时阅读古典老师关于写作这一块的知识点。他提到了10点: 1、你为什么要写作。 2、写作的3种功能。...
    CL生涯规划师阅读 219评论 0 0
  • 刚到福州的第一天下午,在酒店门口巧遇了正面管教联合创始人,鼓励咨询创始人Lynn。记得个把月前,我在wecat上跟...
    吴陈芬Grace阅读 551评论 0 5
  • 声明:文中市场数据来源于网络,主要引用自《车企赞助真人秀 烧钱难挣吆喝》,那篇文章提供了车企赞助的投资与回报的结果...
    在下曾美丽阅读 312评论 0 0