如何自己去写一个鼠标驱动_为什么要用哈密顿采样器（Hamiltonian Monte Carlo），以及如何自己写一个...

背景介绍：（了解采样的可以跳过）

1）为什么需要采样：

简单的分布，比如高斯、exponential、gamma等等的样本都可以直接用numpy.random生成，但复杂的分布需要采样器生成。在贝叶斯、概率编程里面，有很多复杂的分布，而贝叶斯更新需要对这些复杂分布再采样。

2）最简单的3种采样：

CDF采样：（把pdf转成cdf，然后对0-1区间均匀采样）

拒绝采样：

蓝线是真实分布，红线是一个处处y坐标值大于蓝线的函数（正比于一个简单分布，比如均匀分布，或者高斯分布，或者一个分段函数），对红线的分布采样得到样本x，然后计算蓝线(x)/红线(x)作为接受这个样本的概率。

重要性采样

类似拒绝采样，不需要红线函数处处大于蓝线，但每个样本x的权重是蓝线(x)/红线(x)

3）以上三种采样的问题：

cdf采样需要一个数值积分，所以在高维空间采样的时候计算量太大
拒绝采样需要先构造一个总是值大于自己的函数
重要性采样在对分布情况不太清楚的时候采样方差可能很大，即大部分样本的权重都很低
都或多或少需要对分布本身有一定了解，在高维空间不够高效

4）Metropolis-Hastings：

特点：用一个多元正态分布随机游走，不需要对分布本身有太多了解，高维空间时有效，分布本身不需要积分=1，基本只需要保证f(x)在-∞和+∞处=0，0<f(x)<upper bound就行了。

方法：初始点=>随机游走=>如果

则接受这个新的样本，否则按照

的概率接受这个采样。

问题：随机游走的效率仍然不够高。

代码：

def binormal_draw(xprev,beta=1):mean = [0, 0]cov = [[beta**2,0],[0,(1.5*beta)**2]]binormal = scipy.stats.multivariate_normal(mean,cov)return xprev + binormal.rvs()def metropolis(F, qdraw, nsamp, x_init, burnin, thinning=2, beta=1):samples=np.empty((nsamp,2))x_prev = x_initaccepted = 0j = 0for i in range((nsamp+burnin)*thinning):x_star = qdraw(x_prev, beta)logp_star = np.log(F(x_star[0], x_star[1]))logp_prev = np.log(F(x_prev[0], x_prev[1]))logpdfratio_p = logp_star-logp_prevu = np.random.uniform()if np.log(u) <= logpdfratio_p:x_prev = x_starif i >= burnin*thinning and i%thinning == 0:samples[j] = x_starj += 1accepted += 1else:#we always get a sampleif i >= burnin*thinning and i%thinning == 0:samples[j]= x_prevj += 1return samples, accepted

可视化：

http://elevanth.org/blog/2017/11/28/build-a-better-markov-chain/elevanth.org

HMC：

1）直观理解原理：把MCMC的样本随机游走改成在一个势场中具有动能的质点，势场由分布函数f(x)决定，初始动能随机给定，这个质点在运动时间t后的位置记录下来，并且作为下次采样的初始点。

2）相比MH的优点：更新样本点的时候使用了梯度，所以对复杂分布采样比随机试的速度要快。

2）公式、可视化：参考下列文章

Markov Chains: Why Walk When You Can Flow?elevanth.orgXinyu Chen：如何简单地理解「哈密尔顿蒙特卡洛 (HMC)」？zhuanlan.zhihu.com

https://theclevermachine.wordpress.com/2012/11/18/mcmc-hamiltonian-monte-carlo-a-k-a-hybrid-monte-carlo/theclevermachine.wordpress.com

MCMC: Hamiltonian Monte Carlo (a.k.a. Hybrid Monte Carlo)theclevermachine.wordpress.com

https://arxiv.org/pdf/1701.02434.pdfarxiv.orgHamiltonian Monte Carlo explained by Alex Rogozhnikovarogozhnikov.github.io

3）代码：

def HMC(F, u0, n_iter, N_iter, h=0.01): # part 1: 初始化一个orbit的变量，第一行是初始点u0的坐标orbit = torch.zeros((n_iter+1, 2))orbit[0] = u0.detach()u = orbit[0].unsqueeze(0)shape = u.size()# part 2: 循环n_iter次，每次初始化一个v0和u0，然后让leapfrog走N_iter次，每次时长h默认=0.01for k in tqdm(range(n_iter)):v0 = torch.randn(size=u0.shape)u0 = torch.randn(size=u0.shape)*3u, v = leapfrog(F, u0, v0, h, N_iter,shape)u0 = orbit[k]a = float(ratio(F, u0, v0, u, v, shape))r = np.random.rand()if r < a:orbit[k+1] = uelse:orbit[k+1] = u0return orbit #[10:, :]# part 3: leapfrog算法，此外还有velocity verlet等等都可以尝试，但不要用euler
def leapfrog(F, u, v, h, N_iter,shape):v = v - h/2 * grad(F, u)for i in range(N_iter-1):u = u + h * vv = v - h * grad(F, u)u = u + h * vv = v - h/2 * grad(F, u)return u, v# part 4: 对函数F求梯度
def grad(F, u):u = u.detach()if u.requires_grad == False:u = u.requires_grad_()output = -torch.log(F(u))output.backward()ugrad = u.grad.squeeze(0)u = u.squeeze(0)return ugrad# part 5: 细节处理
def unsqueeze(tensor, shape):if len(list(tensor.size())) < len(list(shape)):tensor = tensor.unsqueeze(0)return tensor# part 6: 由于leapfrog计算结果仍然不是绝对的稳定，所以接受概率是min{运行前后总能量的比值,1}
#，大多数情况下都≈1
def ratio(F, u0, v0, u, v, shape):u0 = unsqueeze(u0, shape)u = unsqueeze(u, shape)v0 = unsqueeze(v0, shape)v = unsqueeze(v, shape)w0 = - torch.log(F(u0)) + 0.5*torch.mm(v0, torch.t(v0))w1 = - torch.log(F(u)) + 0.5*torch.mm(v, torch.t(v))return torch.exp(w0-w1)

4）调参细节：

一共有5组主要的超参要调，v0的标准差，u0的位置和标准差，时间步长h，总步长数N_iter，burning。

5）效果：

此处暂时无图。。不过采样的速度总体比Metropolis Hastings要稳定，样本质量较高，使用了梯度所以效率高。

6）进阶：NUTS采样

The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carloarxiv.org