MCMC using Hamiltonian dynamics

Posted mtandhj

tags:

篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了MCMC using Hamiltonian dynamics相关的知识，希望对你有一定的参考价值。

目录

Neal R. M. , MCMC Using Hamiltonian Dynamics[J]. arXiv: Computation, 2011: 139-188.

@article{neal2011mcmc,
title={MCMC Using Hamiltonian Dynamics},
author={Neal, Radford M},
journal={arXiv: Computation},
pages={139--188},
year={2011}}

Hamiltonian Monte Carlo-wiki

算法

先把算法列一下.

Input: 初始值(q), 步长(epsilon)与leapfrog迭代次数(L).

  • 令(q^{(0)} = q);
  • 循环迭代直到停止条件满足(以下第(t)步):
  1. 从标准正态分布中抽取(p), (q = q^{{(t-1)}}), (p^{(t-1)} = p).
  2. [ ag{alg.1} p = p- epsilon abla_q U(q) / 2, ]
  3. 重复 (i=1,2,ldots, L):
    • [ ag{alg.2} q = q + epsilon abla_pH(p), ]
    • 如果(i ot = L):
      [ ag{alg.3} p = p- epsilon abla_q U(q). ]
  4. [ ag{alg.4} p = p- epsilon abla_q U(q) / 2, quad p^* = -p, q^*=q. ]
  5. 计算接受概率
    [ ag{alg.5} alpha = min { 1, exp(-H(q^*,p^*) + H(q^*, p^*)}. ]
  6. 从均匀分布(U(0,1))中抽取(u), 如果(u<alpha), (q^{(t)}=q^*), 否则(q^{(t)}=q^{(t-1)}).

output: ({p^{(t)},q^{(t)}}).

注: 1中的标准正态分布不是唯一的, 但是文中选的便是这个. 4中的(p^*=-p)在实际编写程序的时候可以省去.

符号说明

因为作者从物理方程的角度给出几何解释，所以这里给出的符号一般有俩个含义:
|符号| 概率| 物理|
|--|--|--|
| (q) | 随机变量，服从我们所在意的分布 | 冰球的位置|
| (p)| 用以构造马氏链的额外的变量 | 冰球的动量(mv) |
| (U(q)) | ...与我们所在意的分布有关 | 冰球的势能 |
| (K(p)) | ... | 冰球的动能 |
| (H(q, p)) | 与我们所在意的分布有关 | (H(q, p) = U(q) + K(p)) |

Hamilton方程

物理解释

Hamilton方程((i=1,2,ldots, d), (d)是维度):
[ ag{2.1} frac{dq_i}{dt} = frac{partial H}{partial p_i} ]
[ ag{2.2} frac{dp_i}{dt} = -frac{partial H}{partial q_i} ]

这个东西怎么来的, 大概是因为(H(q, p) = U(q)+K(p)), 如果机械能守恒, 那么随着时间(t)的变化(H(q, p))应该是一个常数, 所以其关于t的导数应该是0.
[ frac{dH}{dt} = sum_i [frac{partial H}{partial p_i} frac{dp_i}{dt}+frac{partial H}{partial q_i}frac{dq_i}{dt}] = ( abla_p H)^T dot{p}+( abla_qH)^T dot{q} = 0, ]
其中(dot{p}=(partial{p_1}/partial t, ldots, partial p_d / partial t)).
而((2.1)), 左边是速度(dot{q} = v), 右边
[ abla_p H = abla_p K = abla_p frac{|p|^2}{2m} = frac{p}{m} = v=dot{q}. ]

不过，估计是先有的(2.2)再有的H吧, 就先这么理解吧. 需要一提的是, (K(p))通常定义为
[ K(p)=p^TM^{-1}p/2, ]
其中(M)是对称正定的, 后面我们可以看到, 这种取法与正态分布相联系起来.
此时:
技术图片

一些性质

可逆 Reversibility

映射(T_s: (q(t), p(t)) mapsto (q(t+s), p(t+s)))是一一的，这个我感觉只能从物理的解释上理解啊, 一个冰球从一个点到另一个点, 现在H确定, 初值确定, 不就相当于整个轨迹都确定了吗, 那从哪到哪自然是一一的, 也就存在逆(T_{-s}), 且只需:
[ frac{dq_i}{dt} = -frac{partial H}{partial p_i}, ]
[ frac{dp_i}{dt} = -(-frac{partial H}{partial q_i}). ]

H的不变性

即
技术图片
当然, 因为我们的算法是离散化的, 所以这个性质只是在(epsilon)比较小的时候近似保持.

保体积 Volume preservation

即假设区域(R={(q,p)})在映射(T_t)的作用下为(R_t={(q(t), p(t)}), 则二者的体积相同, 均为(V).

定义
[ v(t)= int_{R_t} dV = int_{R} det( frac{partial T_t}{partial z}) dV, ]
其中(z = (q, p)). 又
[ z(t) = T_tz = z+t J abla H(z) + mathcal{O}(t^2), ]
其中
[ f:=frac{dz}{dt} = J abla H(z), J = left ( egin{array}{ll} 0_{d imes d} & I_{d imes d} -I_{d imes d} & 0_{d imes d} end{array} ight ). ]
所以
[ frac{partial{T_t}}{partial z} = I + frac{partial f}{partial z}t + mathcal{O}(t^2), ]
又对于任意方阵(A)，
技术图片
所以
[ det (frac{partial{T_t}}{partial z} ) = 1 + mathrm{tr}(partial f / partial z) t + mathcal{O}(t^2), ]
且(mathrm{tr}(partial f / partial z)=mathrm{div} f), 于是
[ frac{d v}{ d t} |_{t=0} = int_{R}mathrm{div} f : d V. ]
又
[ mathrm{div}f = mathrm{div} J abla H(z) = J mathrm{div} abla H(z) = sum_{i=1}^d [frac{partial}{partial q}frac{partial H}{partial p_i}-frac{partial}{partial p}frac{partial H}{partial q_i}] = 0. ]
对于(t=t_0), 我们都可以类似的证明(dv(t_0)/dt=0), 所以(v(t))是常数.

这部分的证明参考自

Here

辛 Symplecticness

技术图片

离散化Hamilton方程 leapfrog方法

下面四幅图, 是(U(q)=q^2/2, K(p)=p^2/2), 起始点为((q, p) = (0, 1)).
技术图片

技术图片

Euler‘s method

如果假定(K(p) = sum_{i=1}^d frac{p_i^2}{2m_i}),
技术图片

Modified Euler‘s method

仅有一个小小的变动,
技术图片

Leapfrog method

技术图片
注意到, 在实际编写程序的时候, 除了第一步和最后一步, 我们可以将(p)的俩个半步合并成一步.
另外从右下角的图可以发现, 因为离散化的缘故, (H(q, p))的值是有偏差的. 但是Leapfrog 方法和 modified Euler方法都是保体积的, 因为每步更新都只改变一个量, 可以验证其雅可比行列式为1.

MCMC

概率与Hamiltonian, 正则(canonical)分布

如何将分布于Hamilton方程联系在一起? 假如, 我们关心的是(q)的分布(P(q)), 则我们构造一个容易采样的分布(P(q)),
[ ag{3.2} P(q, p) = frac{1}{Z}exp(-H(q,p)/T), ]
其中(Z)是规范化的常数, (T)一般取1. 从(3.2)中容易得到(H(q,p)). 事实上此时(q, p)是独立的(这么写是说明直接构造(P(q, p))也是可以的), 则可以分别
[ P(q) = frac{1}{Z_1}exp(-U(q)/T), P(p) = frac{1}{Z_2}exp(-K(p)/T). ]
在贝叶斯统计中, 有
[ U(q) = -log [pi (q) L(D|q)], ]
其中(D)为数据, (L)为似然函数, 与文章中不同, 文章中是(L(q|D)), 应该是笔误.

HMC算法

就是开头提到的算法, 但是其中有一些地方值得思考. (alg.4)我们令(p^*=-p), 这一步在实际中是不起作用的, 既然(K(p)=K(-p))而且在下轮中我们重新采样(p), 我看网上的解释是为了理论, 取反这一部分使得proposal是对称的, 是建议分布(g(p^*, q^*|p, q)=g(p, q| p^*, q^*))? 不是很懂.

关于这个问题的讨论

有点明白了, 首先因为Leapfrog是确定的, 所以(P(q^*, p^*|q, p))非0即1:
[ P(q^*, p^*|q, p) = delta(q^*,p^*,q,p) = left { egin{array}{ll} 1, & T_{Lepsilon} (q, p) = q^*, p^*, , & T_{Lepsilon} (q, p) ot = q^*, p^*. end{array} ight. ]
为了(P(q^*, p^*|q, p)=P(q, p|q^*, p^*)), 如果不取反肯定不行, 因为他就会往下走, 取反的操作实际上就是在可逆性里提到的, 在同样的操作下, (q^*, p^*)会回到(q, p). 于是MH接受概率就退化成了M接受概率. 但是前文也提到了, 取反的操作, 只有在(K(p)=K(-p))的情况下是成立的.

HMC保持正则分布不变的证明 detailed balance

假设({A_k})是((q, p) in mathbb{R}^{2d})空间的一个分割, 其在L次leapfrog的作用, 及取反的操作下的像为({B_k}), 由于可逆性, ({B_k})也是一个分割, 且有相同的体积(V)(保体积性)，则
[ P(A_i) T(B_j|A_i) = P(B_j)T(A_i|B_j). ]
实际上(i ot =j)是时候是显然的, 因为二者都是0. 因为(H)是连续函数, 当(V)变得很小的时候, (H)在(X)区域上的值相当于常数(H_X), 于是
技术图片

所以(3.7)成立.

Detailed balance:

技术图片
其中(R(X))是当前状态属于(X), 拒绝提议的((q^*,p^*))的概率. 注意(sum_{i} T(A_i|B_k)=T(A_k|B_k)=1-R(B_k)). 看上面的连等式可能会有点晕, 注意到, 左端实际上是概率(P{q^{(t)}, p^{(t)}) in B_k}), 最右端是(P{(q^{(t-1)}, p^{(t-1)}) in B_k}), 这样就能明白啥意思了.

遍历性 Ergodicty

马氏链具有遍历性才会收敛到一个唯一的分布上(这部分不了解), HMC是具有这个性质的, 只要(L)和(epsilon)选的足够好. 但是如果选的不过也会导致坏的结果, 比如上面的图, (p^2+q^2=1), 如果我们选择了(Lepsilon approx 2pi), 那么我们的Leapfrog总会带我们回到原点附近, 这就会导致比较差的结果.

HMC的一个例子及优势

下图是:
技术图片
(epsilon=0.25, L=30, q = [-1.5, -1.55]^T, p=[-1, 1^T]).
技术图片

相关系数由(0.95)改为(0.98), (epsilon=0.18, L=20), 随机游走取的协方差矩阵为对角阵, 标准差为(0.18), HMC生成(p)的为标准正态分布.

技术图片

文章中提到, HMC较Randomwalk的优势在于，Randomwalk对协方差很敏感, 而且太大会导致接受率很低, 太小俩俩之间的相关性又会太高.

HMC在实际中的应用和理论

线性变换

有些时候, 我们会对变量施加线性映射(q'=Aq)((A)非奇异方阵), 此时新的密度函数(P'(q') = P(A^{-1}q') / |det (A)|), 其中(P(q))是(q)的密度函数, 相应的我们需要令(U'(q')=U(A^{-1}q')).

如果我们希望线性变化前后不会是的情况变得"更糟", 一个选择是(p'=(A^T)^{-1}p)，则(K'(p')=K(A^Tp')), 如果(K(p) = p^TM^{-1}p / 2), 则
技术图片
其中(M' = (AM^{-1}A^T)^{-1}=(A^{-1})^TMA^{-1}). 此时((q', p'))的更新会使得原先的((q,p))的更新保持, 即
技术图片
所以((q, p))本质上按照原来的轨迹发生着变化.

设想, 我们对(q)的协方差矩阵有一个估计(Sigma), 且近似服从高斯分布, 我们可以对其做Cholesky分解(Sigma=LL^T), 并且令(q'=L^{-1} q), 则(q')的各分量之间就相互独立了, 那么我们很自然的一个选择是(K(p)=p^Tp/2), 那么(q')的各分量的独立性能够保持.

另一个做法是, 保持(q)不变, 但是(K(p) = p^T Sigma p / 2), 此时(q'=L^{-1}q, p'=(L^{T})^{-1}p), 则相当于
[ K(p')=(p')^T{M'}^{-1}p', quad M' = (L^{-1}LL^T(L^{-1})^T)^{-1}=I. ]
所以俩个方法是等价的.

HMC的调整(epsilon, L)

HMC对(epsilon, L)的选择比较严苛.

预先的实验

我们可以对一些(epsilon,L)进行实验, 观察轨迹, 虽然这个做法可能产生误导, 另外在抽样过程中随机选择(epsilon, L)是一个不错的选择.

stepsize (epsilon)

(epsilon)的选择很关键, 如果太大, 会导致低的接受率, 如果太小, 不仅会造成大量的计算成本, 且如果此时(L)也很小, 那么HMC会缺乏足够的探索.

考虑下面的例子:

技术图片
每一次leapfrog将((q(t), p(t)))映射为((q(t+s), p(t+s))), 则
技术图片
((q, p))是否稳定, 关键在于系数矩阵的特征值的模(?还是实部)是否小于1, 特征值为

技术图片
当(epsilon>2sigma)的时候，(4.6)有一个实的大于1的特征值, 当(epsilon < 2 sigma)的时候, 特征值是复制, 且模为:
技术图片
所以, 我们应当选择(epsilon < 2sigma).

在多维问题中, (K(p)=p^Tp/2),如果(q<0)且协方差矩阵为(Sigma), 我们可以取协方差矩阵的最小特征值(非零?)作为步长. 如果(K(p)=p^TM^{-1}p/2), 我们可以通过线性变化将其转换再考虑.

tracjectory length (L)

如何选择(L)也是一个问题, 我们需要足够大的(L)使得每次的探索的足够的, 以便能够模拟出独立的样本, 但是正如前面所讲, 大的(L)不仅会带来计算成本, 而且可能会导致最后结果在起点附近(由周期性带来的麻烦). 而且(L)没法通过轨迹图正确的选择. 一个不错的想法是在一个小的区间内随机选择(L), 这样做可能会减少由于周期性带来的麻烦.

多尺度

我们可以利用(q)的缩放信息, 为不同的(q_i)添加给予不同的(epsilon_i). 比方说在(K(p)=p^Tp/2)的前提下, 应该对(q_i)放大(s_i)倍, 即(q'= q / s_i)((p)不变).

等价的, 可以令(K(p) = p^TM^{-1}p / 2, m_i = 1/ s_i^2)((q)不变), 相当于(q_i'=q_i/s_i, p'=s_ip), 则
(m_i' = s_i (1/ s_i^2)s_i=1), 所以(K(p')=(p')^Tp'/2). 这么做就相当于一次leapfrog为:

技术图片

结合HMC与其它MCMC

当我们所关心的变量是离散的, 或者其对数概率密度((U(q)))的导数难以计算的时候, 结合其它MCMC是有必要的.

Scaling with dimensionality

(U(q) = sum u_i(q_i))的情况

如果(U(q)=sum u_i(q_i))， 且(u_i)之间相互独立(?), 这种假设是可行的, 因为之前已经讨论过， 对于(q)其协方差矩阵为(Sigma), 我们可以通过线性变化使其对角化, 且效能保持.

Cruetz指出, 任何的Metropolis形式的算法在采样密度函数(P(x)=frac{1}{Z}exp(-E(x)))的时候都满足:
技术图片
其中(x)表现在的状态, 而(x^*)表提议. 则根据Jensen不等式可知:
[ 1 = mathbb{E}[exp(Delta)] ge exp(mathbb{E}[-Delta]), ]
所以(mathbb{E}[-Delta]le 0),
[ ag{4.18} mathbb{E}[Delta] ge 0. ]

在(U(q) = sum u_i(q_i))的情况下, 令(Delta_1:=E(x^*)-E(x)), (x=q_i, E(x)=u_i(q_i))或者(x=(q_i, p_i), E(x)=u_i(q_i)=p_i^2/2). 对于整个状态, 我们则用(Delta_d)表示, 则(Delta_d)是所有(Delta_1)的和. 既然(Delta)的平均值均为正, 这会导致接受概率(min (1, exp(-Delta_d))的减小(随着维度的增加), 除非以减小步长作为代价, 或者建议分布的宽度进一步降低(即(x,x^*)尽可能在一个区域内).

因为(exp(-Delta_1) approx 1 -Delta_1+Delta_1^2 / 2), 再根据(4.17)得:
[ ag{4.19} mathbb{E} [Delta_1] approx mathbb{E} [Delta_1^2]/2. ]

故(Delta_1)的方差约是均值的两倍((Delta_1)足够小的时候), 类似的也作用与(Delta_d). 为了有一个比较好的接受率, 我们应当控制(Delta_d)的均值在1左右(小于?), 此时(exp(-1)approx0.3679).

HMC的全局误差(标准差)在(mathcal{O}(epsilon^2))级别, 所以(Delta_1^2)应当在(epsilon^4)级别, 所以(mathbb{E}[Delta_1])也应当在(epsilon^4)级别, 则(mathbb{E}[Delta_d])在(depsilon^4)级别上, 所以为了保持均值为1左右, 我们需要令(epsilon)正比于(d^{1/4}), 相应的(L)为(d^4).

文章中还有关于Randomwalk的分析, 这里不多赘述了.

HMC的扩展和变种

这个不一一讲了, 提一下

分割

假如(H)能够分解成:
技术图片
那么我们可以一个一个来处理, 相当于(T_{1, epsilon} circ T_{2, epsilon}, circ cdots circ T_{k-1, epsilon} circ T_{k, epsilon}). 这个做法依旧是保体积的(既然每一个算子都是保体积的), 但是如果希望其可逆(对称), 这就要求(H_i(q, p)=H_{k-i+1}(q, p).) 这个有很多应用场景:

比如(U(q) = U_0(q)+U_1(q)), 其中计算(U_0)如果是容易的, 我们可以作如下分解:
技术图片
此时有(3M+2)项, (H_1(q, p) = H_k(q, p)=U_1(q)/2), (H_{3i-1}=H_{3i+1}=U_0(q)/2M), (H_{3i}=K(p)/M). 此时对于中间部分, 相当于步长变小了，误差自然会小.

当处理大量数据, 并用到似然函数的时候：
技术图片
不过文章中说这个分解是不对称的, 可明明是对称的啊.

处理约束

有些时候, 我们对(q)有约束条件, 比方(q >v, q<w)等等, 直接给算法:

技术图片

Langevin method, LMC

假设我们使用(K(p)=(1/2)sum p_i^2), (p)从标准正态分布中采样, 每次我们只进行一次leapfrog变计算接受概率, 即
技术图片
以及其接受概率:
技术图片

Windows of states

这个方法试图将(H)的曲线平滑来提高接受率.

我们可以通过任意((q, p))构造一列([(q_0, p_0), ldots, (q_{W-1}, p_{W-1})]). 首先从({0, 1, ldots, W-1})中等概率选一个数, 这个数代表((q, p))在序列中的位置, 记为((q_s, p_s)), 则其前面的可以通过leapfrog ((-epsilon))产生, 后面的通过leapfrog((+epsilon))产生. 所以任意列的概率密度为:

技术图片

然后从((q_{W-1}, p_{W-1}))出发, 通过L-W+1步leapfrog((+epsilon))获得([q_W,p_W,ldots, (q_L, p_L)])并定义提议序列为([(q_L, -p_{L}), ldots, (q_{L-W+1}, p_{L-W+1)}])，计算接受概率:
技术图片
其中(P(q, p)propto exp(-H(q,p))). 设拒绝或者接受后的状态为:([(q_0^+, p_0^+), ldots, (q_{W-1}^+, p_{W-1}^+)]), 依照概率
技术图片
抽取((q_e^+, p_e^+)), 这个就是((q,p))后的下一个状态.

文中说这么做分布不变, 因为:
技术图片
如果没理解错, 前面的部分就是((q_e^+, p_e^+))出现在最开始的序列中的概率, 但是中间的接受概率哪里去了？ 总不能百分百接受吧...

代码

import numpy as np
from collections.abc import Iterator, Callable
from scipy import stats
class Euler:
    """
    Euler方法
    """
    def __init__(self, grad_u, grad_k, epsilon):
        self.grad_u = grad_u
        self.grad_k = grad_k
        self.epsilon = epsilon

    def start(self, p, q, steps):
        trajectory_p = [p]
        trajectory_q = [q]
        for t in range(steps):
            temp = p - self.epsilon * self.grad_u(q)  # 更新一步P
            q = q + self.epsilon * self.grad_k(p)  # 更新一步q
            p = temp
            trajectory_p.append(p)
            trajectory_q.append(q)
        return trajectory_p, trajectory_q

    def modified_start(self, p, q, steps):
        """
        启动！
        :param p:
        :param q:
        :param steps:  步数
        :return: p, q
        """
        trajectory_p = [p]
        trajectory_q = [q]
        for t in range(steps):
            p = p - self.epsilon * self.grad_u(q) #更新一步P
            q = q + self.epsilon * self.grad_k(p) #更新一步q
            trajectory_p.append(p)
            trajectory_q.append(q)
        return trajectory_p, trajectory_q
class Leapfrog:
    """
    Leapfrog 方法
    """
    def __init__(self, grad_u, grad_k, epsilon):
        self.grad_u = grad_u
        self.grad_k = grad_k
        self.epsilon = epsilon
        self.trajectory_q = []
        self.trajectory_p = []

    def start(self, p, q, steps):
        self.trajectory_p.append(p)
        self.trajectory_q.append(q)
        for t in range(steps):
            p = p - self.epsilon * self.grad_u(q) / 2
            q = q + self.epsilon * self.grad_k(p)
            p = p - self.epsilon * self.grad_u(q) / 2
            self.trajectory_q.append(q)
            self.trajectory_p.append(p)
        return p, q

class HMC:
    """
    HMC方法 
    start为进行一次
    accept_prob: 计算接受概率
    hmc: 完整的过程
    acc_rate: 接受概率
    """
    def __init__(self, grad_u, grad_k, hamiltonian, epsilon):
        assert isinstance(grad_u, Callable), "function needed..."
        assert isinstance(grad_k, Callable), "function needed..."
        assert isinstance(hamiltonian, Callable), "function needed..."
        self.grad_u = grad_u
        self.grad_k = grad_k
        self.hamiltonian = hamiltonian
        self.epsilon = epsilon
        self.trajectory_q = []
        self.trajectory_p = []

    def start(self, p, q, steps):
        self.trajectory_p.append(p)
        self.trajectory_q.append(q)
        p = p - self.epsilon * self.grad_u(q) / 2
        for t in range(steps-1):
            q = q + self.epsilon * self.grad_k(p)
            p = p - self.epsilon * self.grad_u(q)
        q = q + self.epsilon * self.grad_k(p)
        p = p - self.epsilon * self.grad_u(q) / 2
        p = -p
        return p, q

    def accept_prob(self, p1, q1, p2, q2):
        """
        :param p1: 原先的
        :param q1:
        :param p2: 建议的
        :param q2:
        :return:
        """
        p1 = np.exp(self.hamiltonian(p1, q1))
        p2 = np.exp(self.hamiltonian(p2, q2))
        alpha = min(1, p1 / p2)
        return alpha

    def hmc(self, generate_p, q, iterations, steps):
        assert isinstance(generate_p, Iterator), "Invalid generate_p"
        self.trajectory_q = [q]
        p = next(generate_p)
        self.trajectory_p = [p]
        count = 0.
        for t in range(iterations):
            tempp, tempq = self.start(p, q, steps)
            if np.random.rand() < self.accept_prob(p, q, tempp, tempq):
                p = tempp
                q = tempq
                self.trajectory_p.append(p)
                self.trajectory_q.append(q)
                count += 1
            p = next(generate_p)
        self.acc_rate = count / iterations
        return self.trajectory


    @property
    def trajectory(self):
        return np.array(self.trajectory_p),                np.array(self.trajectory_q)

class Randomwalk:
    """
    walk: 完整的过程, 实际上Metropolis更新似乎就是start中steps=1, 
    一开始将文章的意思理解错了, 不过将错就错, 这样子也能增加一下灵活性.
    """
    def __init__(self, pdf, sigma):
        assert isinstance(pdf, Callable), "function needed..."
        self.pdf = pdf
        self.sigma = sigma
        self.trajectory = []

    def start(self, q, steps=1):
        for t in range(steps):
            q = stats.multivariate_normal.rvs(mean=q,
                                              cov=self.sigma)
        return q

    def accept_prob(self, q1, q2):
        """
        :param q1: 原始
        :param q2: 建议
        :return:
        """
        p1 = self.pdf(q1)
        p2 = self.pdf(q2)
        alpha = min(1, p2 / p1)
        return alpha
    
    def walk(self, q, iterations, steps=1):
        self.trajectory = [q]
        count = 0.
        for t in range(iterations):
            temp = self.start(q, steps)
            if np.random.rand() < self.accept_prob(q, temp):
                q = temp
                count += 1
            self.trajectory.append(q)
        self.acc_rate = count / iterations
        return np.array(self.trajectory)

以上是关于MCMC using Hamiltonian dynamics的主要内容，如果未能解决你的问题，请参考以下文章

PAT 1122 Hamiltonian Cycle

PAT 1122 Hamiltonian Cycle

PAT 1122 Hamiltonian Cycle

1122 Hamiltonian Cycle (25 分)

PAT 1122 Hamiltonian Cycle[比较一般]

PAT甲级——A1122 Hamiltonian Cycle25