背景介紹：（了解采樣的可以跳過）

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.comhttps://theclevermachine.wordpress.com/2012/11/18/mcmc-hamiltonian-monte-carlo-a-k-a-hybrid-monte-carlo/?theclevermachine.wordpress.comMCMC: Hamiltonian Monte Carlo (a.k.a. Hybrid Monte Carlo)?theclevermachine.wordpress.comhttps://arxiv.org/pdf/1701.02434.pdf?arxiv.orgHamiltonian Monte Carlo explained by Alex Rogozhnikov?arogozhnikov.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 Carlo?arxiv.org

總結

以上是生活随笔為你收集整理的如何自己去写一个鼠标驱动_为什么要用哈密顿采样器（Hamiltonian Monte Carlo），以及如何自己写一个...的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。