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-ρ²)^-n/2 exp{(2*(1-ρ²))) ^-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（θ₀,θ₁…θ_t-1）+sdεId，t>=N₀}
C_t是后面不断变化的协方差矩阵
s_d是一个比例因子，仅与参数n有关，以保证α在合适范围。
sd=2.4²/n;cov(θ₀,θ₁…θ_t-1)是历史样本协方差
cov(θ₀,θ₁…θ_t-1)=（t-1）^-1 (sum(θ_iθ_i^T )¹ -t * θ_t-1θ_t-1^T )
ε是一个较小的正数，以保证C_t的非奇异，一般取10^-5 1

7 DR算法

该算法是对建议分布进行局部自适应，以提高采样效率
其基本思想：为了提高新样本接受概率，当u大于计算接受概率α₁时，不立即拒绝候选样本，而是进入第二层建议候选点并计算接受概率α₂，进行在判断
达到目的：候选样本代替新样本的可能性大大提高
减少采样中经常出现的样本平滑段情况
具体过程：假设马尔科夫链的当前位置θ_n=θ，依据标准MH从建议分布
q₁(θ,x₁)选取一个候选点x₁,则此次的接受概率α₁：

一旦此候选点被拒绝，为不使θ_n+1=θ发生，因此又依据第二层的建议分布
q₂（θ，x₁，.），不仅取决于当前链的位置，而且也与被拒绝的候选点x₁有关，则第二层的接受概率α₂：
第二层接受概率α₂（θ，x₁ ,x₂）=

如此循环可延缓拒绝至第i层，但会比较复杂，超过2就不好编了。第3层勉强可以。假设第i层建议分布为q_i（θ，x₁，…x_i-1）,则其接受概率：

至此，所有接受概率均已计算出，是目标分布π的可逆性在每一层延缓中被独立保存下来，延缓拒绝的过程可以在任一层被打断，因此我们可以事先设定一个最大延缓数，比如5层了会被拒绝，就接受链的之前位置。
不难看出Dr只是针对某一次抽样，使其走向被接受，以免被浪费，是一种局部调整的思想。

8 DRAM算法

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

9 参考文献

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

---

1. i从0到t-1 ↩︎