MCMC抽样学习记录(一)

1.伪随机数

首先你要弄懂伪随机数如何生成
Inside_Zhang
一般是线性的 还有高斯的
计算方法有 Box–Muller_transform (an efficient polar form),Ziggurat algorithm 等
验证的话可以计算这个序列的 entropy 或者用压缩算法计算该序列的冗余。
求序列的熵 程序实现
熵的概念以及理论求法
近似熵 样本熵 模糊熵
时间序列

2 各种采样方法

逆采样,重要性采样
概率值融合 把对一个未知分布的函数转到一个已知简单分布(均匀分布)的简单函数。LOTUS
拒绝抽样 : 高维的情况下,Rejection Sampling 两个问题,寻找合适的 q 分布和 k 值。这两个问题会导致拒绝率很高,无用计算增加。

3马尔科夫链

自己定一个法则,比如最常用的例子 上中下层阶级的跃迁问题如果下中上分别对应1,2,3
满足这个矩阵[x11=0.65 x12=0.28 x13=0.07 ;x21=0.15 x22=0.67 x23=0.18 ;x31=0.12 x32=0.36 x33=0.52]
假设开始都分布在下层[1 0 0]
经过11代左右可以得到比较平稳的分布[0.2870 0.4883 0.2247]
马氏链的收敛定理十分重要

4 gibbs采样

Gibbs Sampling 就是从条件概率中选择一个变量(分子),然后对该变量(分子)进行采样。gibbs 采样是已知权重参数和一个 v 变量,通过采样得到 h。通过 h 采样又可以得到另一个 v ,如此交替采样,就可以逐渐收敛于联合分布了

5 M-H算法

M-H算法的效率与建议分布q的选择有关。他的影响有两个方面,一个是接受率,而另一个是由链所覆盖的样本空间的区域——Chib and Greenberg
另一方面,如果扩展域太大,生成的备选值的接受概率会很低,反之则链需要较长时间遍历密度的支撑。关于q1(随机游动建议密度)的结果表明,若目标分布与建议密度都为正态,(1)一维问题接受率取0.45(2)维数接近无穷取0.23(3)维数为6时最佳接受概率为0.25
对于独立链,取q(θ,θ‘)=q(θ’)最重要的是保证建议密度q(θ’)的尾部可以控制目标函数的尾部
协方差矩阵可以用来表示多维随机变量的概率密度
设D={yi=(y1i,y2i)’,i=1,2,…n}是从二元正态分布N2(0,∑)
∑=(1 ρ
ρ 1)
ρ后验密度
π(ρ|D)∝(1-ρ2)-n/2 exp{(2*(1-ρ2))) -1 (S11-2ρS12+S22)
其中ρ∈(-1,1),Srs=∑yri+ysi, r,s=1,2

在标准MH算法,终止运算的条件是:1经过若干次迭代,围绕随机变量θ=θ的估计值小幅波动;2 被估计参数的均值趋于稳定;3 直接设置迭代次数
对于修正后的有限元模型,通常采用一定的准则来评价是否真的反映实际结构的响应,常用的准则包括频率相关性分析、频率与振型相关性分析频响函数相关分析和交叉正交性相关分析等。

6 AM算法

adapt metropolis 自适应metropolis
全局自适应
建议分布很关键 理论上可以是任何形式,但公认的事实是建议分布和目标分布符合程度越高,模拟的效果越好,采样效率随之提高。
选取过大或过小的协方差矩阵都不能得到理想的采样效果,需要寻找协方差矩阵的方法,自主选择建议分布
首先依据先验知识构建初始高斯建议分布,然后利用已产生的markov链的全部历史样本来调整高斯分布的协方差矩阵,从而改变建议分布的空间延展性和方向性,自适应地逼近目标分布

具体做法:
假定非适应阶段的长度为N0,在此阶段中高斯建议分布的协方差矩阵依据经验取为固定值C0
Ct{C0,t<N0;sd * cov(θ01…θt-1)+sdεId,t>=N0}
Ct是后面不断变化的协方差矩阵
sd是一个比例因子,仅与参数n有关,以保证α在合适范围。
sd=2.42/n;cov(θ01…θt-1)是历史样本协方差
cov(θ01…θt-1)=(t-1)-1 (sum(θiθiT )1 -t * θt-1θt-1T )
ε是一个较小的正数,以保证Ct的非奇异,一般取10-5 1
1

7 DR算法

该算法是对建议分布进行局部自适应,以提高采样效率
其基本思想:为了提高新样本接受概率,当u大于计算接受概率α1时,不立即拒绝候选样本,而是进入第二层建议候选点并计算接受概率α2,进行在判断
达到目的:候选样本代替新样本的可能性大大提高
减少采样中经常出现的样本平滑段情况
具体过程:假设马尔科夫链的当前位置θn=θ,依据标准MH从建议分布
q1(θ,x1)选取一个候选点x1,则此次的接受概率α1
3.4
一旦此候选点被拒绝,为不使θn+1=θ发生,因此又依据第二层的建议分布
q2(θ,x1,.),不仅取决于当前链的位置,而且也与被拒绝的候选点x1有关,则第二层的接受概率α2
第二层接受概率α2θ,x1 ,x2)=
在这里插入图片描述
如此循环可延缓拒绝至第i层,但会比较复杂,超过2就不好编了。第3层勉强可以。假设第i层建议分布为qi(θ,x1,…xi-1),则其接受概率:
在这里插入图片描述
至此,所有接受概率均已计算出,是目标分布π的可逆性在每一层延缓中被独立保存下来,延缓拒绝的过程可以在任一层被打断,因此我们可以事先设定一个最大延缓数,比如5层了会被拒绝,就接受链的之前位置。
不难看出Dr只是针对某一次抽样,使其走向被接受,以免被浪费,是一种局部调整的思想。

8 DRAM算法

delay rejecting adaptive metropolis method
6、7两者结合而成,既考虑全局也考虑局部。
随机游走DRAM算法的具体步骤
(1)假设n时刻,markov链的当前位置为θn,依据第一层建议分布产生候选点x1=N(·|θn,covn);
(2)根据式3.4计算接受概率α1
(3)在[0,1]随机产生一个u,若α1>=u则接受候选点θn+1=x1,→ (5),否则→ (4)
(4)依据第i层的建议分布产生候选点xi=Ni(·|θn,x1,x2,…,xi-1,covni);
(5)进入下一时刻n+1,若n+1>=N0,则进入自适应阶段,按照式 (3.1) 自主调整协方差covn+1,否则就取covn+1=C0
(6)重复以上步骤,直至产生足够的样本。

9 参考文献

[1]基于贝叶斯模型更新的结构损伤识别方法改进及应用.张建新.重庆大学.


  1. i从0到t-1 ↩︎

已标记关键词 清除标记
相关推荐
©️2020 CSDN 皮肤主题: 大白 设计师:CSDN官方博客 返回首页