采样-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)})\)
