[MCSM] Slice Sampler
1. 引言
之前介绍的MCMC算法都具有一般性和通用性(这里指Metropolis-Hasting 算法),但也存在一些特殊的依赖于仿真分布特征的MCMC方法。在介绍这一类算法(指Gibbs sampling)之前,本节将介绍一种特殊的MCMC算法。 我们重新考虑了仿真的理论基础,建立了Slice Sampler。
考虑到[MCSM]伪随机数和伪随机数生成器中提到的产生服从f(x)密度分布随机数等价于在子图f上产生均匀分布,即
类似笔记“[MCSM] Metropolis-Hastings 算法”(文章还没写好),考虑采用马尔可夫链的稳态分布来等价上的均匀分布,以此作为f分布的近似。很自然的想法是采用随机行走(random walk)。这样得到的稳态分布是在集合上的均匀分布。
2. 2D slice sample
有很多方法实现在集合上的"random walk",最简单的就是一次改变一个方向上的取值,每个方向的改变交替进行,由此得到的算法是 2D slice sampler
在第t次迭代中,执行
1.
2. , 其中
举例
其中,是归一化因子,代码如下,第一幅图是前10个点的变化轨迹,第二幅图表明初始点的选取影响不大
% p324 T = 0:10000; T = T/10000; % N(3,1) y = exp(-(T+3).^2/2); plot(T,y); hold on; x = 0.25; u = rand *(exp(-(x+3).^2/2)); x_s = [x]; u_s = [u]; for k = 1:10;limit = -3 + sqrt(-2*log(u));limit = min([limit 1]);x = rand * limit;x_s = [x_s x];u_s = [u_s u];u = rand *(exp(-(x+3).^2/2));x_s = [x_s x];u_s = [u_s u]; end plot(x_s,u_s,'-*'); hold off;%% x = 0.01; u = 0.01; x_s = [x]; u_s = [u]; for k = 1:50;limit = -3 + sqrt(-2*log(u));limit = min([limit 1]);x = rand * limit;x_s = [x_s x];u_s = [u_s u];u = rand *(exp(-(x+3).^2/2));x_s = [x_s x];u_s = [u_s u]; end figure; subplot(1,3,1); plot(x_s,u_s,'*');hold on;plot(T,y); x = 0.99; u = 0.0001; x_s = [x]; u_s = [u]; for k = 1:50;limit = -3 + sqrt(-2*log(u));limit = min([limit 1]);x = rand * limit;x_s = [x_s x];u_s = [u_s u];u = rand *(exp(-(x+3).^2/2));x_s = [x_s x];u_s = [u_s u]; end subplot(1,3,2); plot(x_s,u_s,'*');hold on;plot(T,y); x = 0.25; u = 0.0025; x_s = [x]; u_s = [u]; for k = 1:50;limit = -3 + sqrt(-2*log(u));limit = min([limit 1]);x = rand * limit;x_s = [x_s x];u_s = [u_s u];u = rand *(exp(-(x+3).^2/2));x_s = [x_s x];u_s = [u_s u]; end subplot(1,3,3); plot(x_s,u_s,'*');hold on;plot(T,y);
View Code
3. General Slice Sampler
有时候面临的概率密度函数不会那么简单,此时面临的困难主要在于无法在第二次更新的时候找到集合的范围。但有时我们可以将概率密度函数分解为多个较为简单的函数之积,即
Slice Sampler
1.
2. ,其中
看着挺高级好用的,实际上也只是能用的,一是本身就很难求,第二是即使求出来了,这个满足均匀分布的变量也很难得到,比如说书上的例子(Example 8.3)
很自然的分成了
但是集合完全没有办法用,求其中一个,然后拒绝不满足要求的看起来是可行的,但是效率实在是太低了(效率低实际上是我写错了,实际上还可以)
代码如下(代码是MATLAB的,画出来的图不好看,这个图是作者的R代码画出来的)
x = 0; u1 = rand*(1+sin(3*x)^2); u2 = rand*(1+cos(5*x)^4); u3 = rand*(exp(-x^2/2)); x_s = zeros(1,10000); for k = 1:10000limit = sqrt(-2*log(u3));x = -limit + 2*limit*rand;while((sin(3*x))^2<u1-1 || (cos(5*x))^4<u2-1)x = -limit + 2*limit*rand;endu1 = rand*(1+sin(3*x)^2);u2 = rand*(1+cos(5*x)^4);u3 = rand*(exp(-x^2/2));x_s(k) = x; end hist(x_s,100);
View Code
4. 收敛性
不会
转载于:https://www.cnblogs.com/sea-wind/p/4531623.html
[MCSM] Slice Sampler相关推荐
- 【ML】Markov Chain Monte Carlo(MCMC)---Slice sampler(切片采样)和Hierarchical Models(层次模型)
导航 Slice sampler 2D slice sample General Slice Sampler Hierarchical models python Code download Refe ...
- 贝叶斯集锦:从MC、MC到MCMC
转载自: #####一份草稿 贝叶斯计算基础 一.从MC.MC到MCMC 斯坦福统计学教授Persi Diaconis是一位传奇式的人物.Diaconis14岁就成了一名魔术师,为了看懂数学家Fell ...
- Markov Chain Monte Carlo
转载至https://zhuanlan.zhihu.com/p/25610149 [数据分析] Markov Chain Monte Carlo Markov Chain Monte Carlo简称M ...
- MCMC: Metropolis-Hastings, Gibbs and slice sampling
从标题中即已看出,Metropolis-Hastings.Gibbs and slice sampling是MCMC(基于马尔科夫链的蒙特卡洛采样算法)的变体及改进算法. With MCMC, we ...
- Rocksdb Slice使用中的一个小坑
本文记录一下使用Rocksdb Slice过程中的一个小小坑,差点没一口老血吐出来. rocksdb的Slice 数据结构是一个小型得不可变类string数据结构,设计出来的目的是为了保证rocksd ...
- 【Go】Go基础(六):数组、切片slice、映射map
1.数组 1.1 数组定义 数组:数组是具有相同 唯一类型 的一组已编号且长度固定的数据项序列. 数组长度必须是一个常量表达式,并且必须是一个非负整数.数组长度也是数组类型的一部分,所以[5]int和 ...
- 【GoLang】深入理解slice len cap什么算法? 参数传递有啥蹊跷?
先上结论 1.内置append函数在现有数组的长度 < 1024 时 cap 增长是翻倍的,再往上的增长率则是 1.25,至于为何后面会说. 2.Go语言中channel,slice,map这三 ...
- Jmeter调试工具---Debug Sampler
使用场景:脚本开发是,调试用(正式测试时需删除),Debug Sampler会把我们自定义的变量输出在response data中 使用设置:JMeter properties和System prop ...
- static slice是什么呢?
参考:Computing Static Slice for Java Programs - 百度学术 Slicing is an analysis technique that reduces pro ...
最新文章
- python tqdm_Python基础 | 一个被忽视的神器tqdm
- SAP LT Replication Server与SAP HANA中与Replication相关的表
- element el-input 自动获取焦点和IE下光标位置解决方法
- 区分Debug版还是Relase版
- 排序中减治法算法伪代码_算法浅谈——分治算法与归并、快速排序(附代码和动图演示)...
- (23)FPGA面试题常用逻辑电平
- 一加到1亿。C语言_一加官方道歉!这下良心了:老用户欢呼
- 网络负载平衡(Network Load Balancing)的工作原理
- 2017CCPC哈尔滨 D:X-Men
- 【比赛】NOIP2017 列队
- linux代码诊断有没有link,Linux下判断网线是否插入的代码
- css 文字溢出文本时省略号代替
- button组件 untiy_Unity 3D Button控件
- 项目保密协议书(范本)
- 动态规划之挖金矿问题(Python and Java)
- matlab画组合立方体,matlab小程序 画立方体
- 计算机通识必修课程学什么内容,计算机通识课程教学平台研究与探索.doc
- Netlink的简介及使用方法
- 集五福主题的微信图文排版攻略已到!
- maya藤蔓插件_3DMax藤蔓生长插件Guruware Ivy For 3dsmax中英文版本