采样-MH采样

https://www.bilibili.com/video/BV17D4y1o7J2
wiki
https://en.wikipedia.org/wiki/Metropolis–Hastings_algorithm
知乎-介绍了为什么pdf(loc=x_start)的好处
https://www.zhihu.com/question/63305712

1.马尔科夫链MCMC采样

2.问题

有了MCMC采样,但是问题是想要得到对应收敛分布的转移矩阵P如何求得?

2.1 细致平稳条件

$\pi(i)P(i,j)=\pi(j)P(j,i)$
$\pi(x)是指系统在i状态的分布函数,个人任务写成\pi_i(x)更合适,这个老师写的不合理,因为x表示随机变量X的取值,比如系统有A,B,C三种取值,假设是0.5,0.3,0.2,\pi_i("A")=0.5表示i状态时X="A"的概率,\pi_i(x)则表示为一个向量,为[0.5,0.3,0.2]$
$P(i,j)是指系统从i状态转移到j状态的状态转移矩阵$

$回到方程本身\pi_i(x)P(i,j)=\pi_j(x)P(j,i),其代表的是系统从i状态转移到j状态 和 j转移到i状态的概率是一样的$
$对应的如果\pi_i数值较高,则对应的P(i,j)概率较低,并且对应的 \pi_j(x) 较低 P(j,i)较高$

2.2 引入参数$\alpha$

$假设我们现在引入任意随机矩阵Q,显然$
$\pi_i(x)Q(i,j) \ne \pi_j(x)Q(j,i)$
$所以我们引入参数\alpha,令其强行相等$
$\pi_i(x)Q(i,j)\alpha(i,j)=\pi_j(x)Q(j,i)\alpha(j,i)$
令
$\alpha(i,j)=\pi_j(x)Q(j,i)\in[0,1],\alpha(i,j)代表的是从i状态转义到j状态的概率(随着状态转义矩阵P)$
$\alpha(j,i)=\pi_i(x)Q(i,j)\in[0,1]$
两边即相等
$我们得到,P(i,j)=Q(i,j)\alpha(i,j)$
$代表Q不等于P,但是却有一定的概率=P$
$所以结论就是$
$1.给定任意转移矩阵Q,然后计算\alpha$
$2.以Q作为转移矩阵采样,并生成随机数[0,1],若随机数<\alpha,则接受采样,否则丢弃继续采样$

2.3 MCMC算法

$具体算法$
$输入任意给定的马尔科夫链状态转移矩阵 Q，目标平稳分布\pi(x) ，设定状态转移次数阈值n_1，需要的样本数n_2$
$从任意简单概率分布得到初始状态值x_0；$
$for\ t =0\ in\ \ n_1+n_2 -1$
$a. 从条件概率分布 Q(x|x_1)得到样本值 x$
$b. 从均匀分布中采样U \sim [0,1]$
$c. 如果 u < \alpha(x_1,x^*)=\pi(x^*)Q(x^*,x_t)，则接受x_t \rightarrow x^* ，即x_{t+1}=x$
$d. 否则不接受转移， t=max\{t-1,0\}$

3.MH算法

在前面的MCMC算法中若$\alpha$过小容易造成采样效率低下,为了解决这个问题MH算法做了优化
$优化的公式也很简单,令\alpha=min\{\frac{\pi_j(x)Q(j,i)}{\pi_i(x)Q(i,j)},1\}$
$采样算法和上面的MCMC一致，只是取代了\alpha的计算公式,证明略$

4.python MH算法实现

点击查看代码

from scipy.stats import norm

def norm_dist_prob(theta):
    y = norm.pdf(theta, loc=3, scale=2)
    return y

T = 5000
pi = [0 for i in range(T)]
sigma = 1
t = 0
while t < T-1:
    t = t + 1
    pi_star = norm.rvs(loc=pi[t - 1], scale=sigma, size=1, random_state=None)   #状态转移进行随机抽样
    alpha = min(1, (norm_dist_prob(pi_star[0]) / norm_dist_prob(pi[t - 1])))   #alpha值
    u = random.uniform(0, 1)
    if u < alpha:
        pi[t] = pi_star[0]
    else:
        pi[t] = pi[t - 1]


plt.scatter(pi, norm.pdf(pi, loc=3, scale=2),label='Target Distribution')
num_bins = 50
plt.hist(pi, num_bins, normed=1, facecolor='red', alpha=0.7,label='Samples Distribution')
plt.legend()
plt.show()

5.Metropolis 算法

上面的算法其实很难理解，究其原因
1.之前的demo是介绍了离散值情况下的MC采样,这个很好理解,但如果是连续值如何写转移矩阵?
2.norm.rvs(loc=pi[t - 1], scale=sigma, size=1, random_state=None) 代码看不懂,关键是算法没理解
3.代码中没看到乘以矩阵Q的操作,没看懂

对于第一点
$对于离散值,矩阵中每个值代表的是Q(X=|Y=)是一个带条件的概率质量函数(pmf)$
$那么扩展到连续纸,转移概率矩阵就是一个带条件的概率密度函数(pdf)$
$所以总结下,若是离散值,那么就写成状态转义矩阵Q(x|y),若是连续值,就写成g(x|y),代表概率密度函数$
$\color{red}{这里注意下一般的书上或者教材上都会写成条件概率Q{x|y},这种数学上表达式没错,但不用关心这个条件y,你也没法求解,脱离这个条件y,他的本质还是某个概率密度函数pdf/pmf,只要你能找到这种你想要的pdf/pmf即可}$
再看第二,三点
先看第一个wiki上的介绍
https://en.wikipedia.org/wiki/Random_walk
$为了取到合适的Q/g,Metropolis算法要求必须满足Q(x|y)=Q(y|x),或者g(x|y)=g(y|x),也就是这个分布是对称的,然后重点来了$
wiki原文

---

A usual choice is to let ${\displaystyle g(x\mid y)}{\displaystyle g(x\mid y)}$be a Gaussian distribution centered at ${\displaystyle y}$, so that points closer to ${\displaystyle y}$ are more likely to be visited next, making the sequence of samples into a random walk [a]. The function ${\displaystyle g}$ is referred to as the proposal density or jumping distribution.

---

再抄一段知乎上的(有版权~~~大家自己去看下)
https://www.zhihu.com/question/63305712
从 第三：意外之喜。 这里开始看
$大致意思就是通过wiki的方法得到了Metropolis算法的精髓$
$若MH算法中的提议分布(就是Q/g)如果是对称的,那么有Q(x|y)=Q(y|x),或者g(x|y)=g(y|x),那么MH算法可以简化为$
$\alpha=min(1,\frac{\pi_j(x)}{\pi_i(x)})$