网站开发能封装成app吗/实体店铺引流推广方法

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$, 步长$ϵ$与leapfrog迭代次数$L$.

• 令$q^{(0)}=q$;
• 循环迭代直到停止条件满足(以下第$t$步):

1. 从标准正态分布中抽取$p$, $q=q^{(t-1)}$, $p^{(t-1)}=p$.
2. $$p=p-ϵ{\nabla }_{q}U(q)/2,\text{\hspace{1em}(alg.1)}$$
3. 重复 $i=1,2,\dots ,L$:
   • $$q=q+ϵ{\nabla }_{p}H(p),\text{\hspace{1em}(alg.2)}$$
   • 如果$i\ne L$:
     $$p=p-ϵ{\nabla }_{q}U(q).\text{\hspace{1em}(alg.3)}$$
4. $$p=p-ϵ{\nabla }_{q}U(q)/2,p^{\ast }=-p,q^{\ast }=q.\text{\hspace{1em}(alg.4)}$$
5. 计算接受概率
   α=min⁡{1,exp⁡(−H(q∗,p∗)+H(q∗,p∗)}.(alg.5)\tag{alg.5} \alpha = \min \{ 1, \exp(-H(q^*,p^*) + H(q^*, p^*)\}. α=min{1,exp(−H(q∗,p∗)+H(q∗,p∗)}.(alg.5)
6. 从均匀分布$U(0,1)$中抽取$u$, 如果$u<\alpha$, $q^{(t)}=q^{\ast }$, 否则$q^{(t)}=q^{(t-1)}$.

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

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

符号说明

因为作者从物理方程的角度给出几何解释，所以这里给出的符号一般有俩个含义:

Hamilton方程

物理解释

Hamilton方程($i=1,2,\dots ,d$, $d$是维度):
$$\frac{d{q}_i}{dt}=\frac{\partial H}{\partial p_i}\text{\hspace{1em}(2.1)}$$
$$\frac{d{p}_i}{dt}=-\frac{\partial H}{\partial q_i}\text{\hspace{1em}(2.2)}$$

这个东西怎么来的, 大概是因为$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{d{p}_i}{dt}+\frac{\partial H}{\partial q_i}\frac{d{q}_i}{dt}]=({\nabla }_{p}H{)}^T\dot{p}+({\nabla }_{q}H{)}^T\dot{q}=0,$$
其中$\dot{p}=(\partial p_1/\partial t,\dots ,\partial p_d/\partial t)$.
而$(2.1)$, 左边是速度$\dot{q}=v$, 右边
$${\nabla }_{p}H={\nabla }_{p}K={\nabla }_p\frac{|p{|}^2}{2m}=\frac{p}{m}=v=\dot{q}.$$

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

一些性质

可逆 Reversibility

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

H的不变性

即

当然, 因为我们的算法是离散化的, 所以这个性质只是在$ϵ$比较小的时候近似保持.

保体积 Volume preservation

即假设区域R={(q,p)}R=\{(q,p)\}R={(q,p)}在映射$T_t$的作用下为Rt={(q(t),p(t)}R_t=\{(q(t), p(t)\}Rt​={(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_{t}z=z+tJ\nabla H(z)+\mathcal{O}(t^2),$$
其中
$$\begin{gathered} f:=\frac{dz}{dt}=J\nabla H(z), \\ J=\left(\begin{array}{ll}{0}_{d\times d}& I_{d\times d}\\ -I_{d\times d}& 0_{d\times d}\end{array}\right). \end{gathered}$$
所以
$$\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{t}\mathrm{r}(\partial f/\partial z)t+\mathcal{O}(t^2),$$
且$\mathrm{t}\mathrm{r}(\partial f/\partial z)=\mathrm{d}\mathrm{i}\mathrm{v}f$, 于是
$$\frac{dv}{dt}{|}_{t=0}={\int }_R\mathrm{d}\mathrm{i}\mathrm{v}fdV.$$
又
$$\mathrm{d}\mathrm{i}\mathrm{v}f=\mathrm{d}\mathrm{i}\mathrm{v}J\nabla H(z)=J\mathrm{d}\mathrm{i}\mathrm{v}\nabla 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)$,
$$P(q,p)=\frac{1}{Z}\exp (-H(q,p)/T),\text{\hspace{1em}(3.2)}$$
其中$Z$是规范化的常数, $T$一般取1. 从(3.2)中容易得到$H(q,p)$. 事实上此时$q,p$是独立的(这么写是说明直接构造$P(q,p)$也是可以的), 则可以分别
$$\begin{gathered} P(q)=\frac{1}{Z_1}\exp (-U(q)/T), \\ P(p)=\frac{1}{Z_2}\exp (-K(p)/T). \end{gathered}$$
在贝叶斯统计中, 有
$$U(q)=-\log [\pi (q)L(D|q)],$$
其中$D$为数据, $L$为似然函数, 与文章中不同, 文章中是$L(q|D)$, 应该是笔误.

HMC算法

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

关于这个问题的讨论

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

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

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

所以(3.7)成立.

Detailed balance:

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

遍历性 Ergodicty

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

HMC的一个例子及优势

下图是:

$ϵ=0.25,L=30,q=[-1.5,-1.55{]}^T,p=[-1,1^T]$.

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

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

HMC在实际中的应用和理论

线性变换

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

如果我们希望线性变化前后不会是的情况变得"更糟", 一个选择是$p^{\prime }=(A^T{)}^{-1}p$，则$K^{\prime }(p^{\prime })=K(A^T{p}^{\prime })$, 如果$K(p)=p^T{M}^{-1}p/2$, 则

其中$M^{\prime }=(A{M}^{-1}{A}^T{)}^{-1}=(A^{-1}{)}^{T}M{A}^{-1}$. 此时$(q^{\prime },p^{\prime })$的更新会使得原先的$(q,p)$的更新保持, 即

所以$(q,p)$本质上按照原来的轨迹发生着变化.

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

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

HMC的调整$ϵ,L$

HMC对$ϵ,L$的选择比较严苛.

预先的实验

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

stepsize $ϵ$

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

考虑下面的例子:

每一次leapfrog将$(q(t),p(t))$映射为$(q(t+s),p(t+s))$, 则

$(q,p)$是否稳定, 关键在于系数矩阵的特征值的模(?还是实部)是否小于1, 特征值为

当$ϵ>2\sigma$的时候，(4.6)有一个实的大于1的特征值, 当$ϵ<2\sigma$的时候, 特征值是复制, 且模为:

所以, 我们应当选择$ϵ<2\sigma$.

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

tracjectory length $L$

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

多尺度

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

等价的, 可以令$K(p)=p^T{M}^{-1}p/2,m_i=1/s_i^2$($q$不变), 相当于$q_i^{\prime }=q_i/s_i,p^{\prime }=s_{i}p$, 则
$m_i^{\prime }=s_i(1/s_i^2)s_i=1$, 所以$K(p^{\prime })=(p^{\prime }{)}^T{p}^{\prime }/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^{\ast }$表提议. 则根据Jensen不等式可知:
$$1=\mathbb{E}[\exp (\Delta )]\ge \exp (\mathbb{E}[-\Delta ]),$$
所以$\mathbb{E}[-\Delta ]\le 0$,
$$\mathbb{E}[\Delta ]\ge 0.\text{\hspace{1em}(4.18)}$$

在$U(q)=\sum u_i(q_i)$的情况下, 令${\Delta }_1:=E(x^{\ast })-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^{\ast }$尽可能在一个区域内).

因为$\exp (-{\Delta }_1)\approx 1-{\Delta }_1+{\Delta }_1^2/2$, 再根据(4.17)得:
$$\mathbb{E}[{\Delta }_1]\approx \mathbb{E}[{\Delta }_1^2]/2.\text{\hspace{1em}(4.19)}$$

故${\Delta }_1$的方差约是均值的两倍(${\Delta }_1$足够小的时候), 类似的也作用与${\Delta }_d$. 为了有一个比较好的接受率, 我们应当控制${\Delta }_d$的均值在1左右(小于?), 此时$\exp (-1)\approx 0.3679$.

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

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

HMC的扩展和变种

这个不一一讲了, 提一下

分割

假如$H$能够分解成:

那么我们可以一个一个来处理, 相当于$T_{1,ϵ}\circ T_{2,ϵ},\circ \cdots \circ T_{k-1,ϵ}\circ T_{k,ϵ}$. 这个做法依旧是保体积的(既然每一个算子都是保体积的), 但是如果希望其可逆(对称), 这就要求$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),\dots ,(q_{W-1},p_{W-1})]$. 首先从{0,1,…,W−1}\{0, 1, \ldots, W-1\}{0,1,…,W−1}中等概率选一个数, 这个数代表$(q,p)$在序列中的位置, 记为$(q_s,p_s)$, 则其前面的可以通过leapfrog ($-ϵ$)产生, 后面的通过leapfrog($+ϵ$)产生. 所以任意列的概率密度为:

然后从$(q_{W-1},p_{W-1})$出发, 通过L-W+1步leapfrog($+ϵ$)获得$[q_W,p_W,\dots ,(q_L,p_L)]$并定义提议序列为$[(q_L,-p_L),\dots ,(q_{L-W+1},p_{L-W+1)}]$，计算接受概率:

其中$P(q,p)\propto \exp (-H(q,p))$. 设拒绝或者接受后的状态为:$[(q_0^{+},p_0^{+}),\dots ,(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_uself.grad_k = grad_kself.epsilon = epsilondef start(self, p, q, steps):trajectory_p = [p]trajectory_q = [q]for t in range(steps):temp = p - self.epsilon * self.grad_u(q)  # 更新一步Pq = q + self.epsilon * self.grad_k(p)  # 更新一步qp = temptrajectory_p.append(p)trajectory_q.append(q)return trajectory_p, trajectory_qdef 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) #更新一步Pq = q + self.epsilon * self.grad_k(p) #更新一步qtrajectory_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_uself.grad_k = grad_kself.epsilon = epsilonself.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) / 2q = q + self.epsilon * self.grad_k(p)p = p - self.epsilon * self.grad_u(q) / 2self.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_uself.grad_k = grad_kself.hamiltonian = hamiltonianself.epsilon = epsilonself.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) / 2for 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) / 2p = -preturn p, qdef 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 alphadef 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 = temppq = tempqself.trajectory_p.append(p)self.trajectory_q.append(q)count += 1p = next(generate_p)self.acc_rate = count / iterationsreturn self.trajectory@propertydef 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 = pdfself.sigma = sigmaself.trajectory = []def start(self, q, steps=1):for t in range(steps):q = stats.multivariate_normal.rvs(mean=q,cov=self.sigma)return qdef accept_prob(self, q1, q2):""":param q1: 原始:param q2: 建议:return:"""p1 = self.pdf(q1)p2 = self.pdf(q2)alpha = min(1, p2 / p1)return alphadef 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 = tempcount += 1self.trajectory.append(q)self.acc_rate = count / iterationsreturn np.array(self.trajectory)