开这个坑慢慢写吧…

Langevin dynamic 和 Hamiltonian Monte Carlo

众所周知，sampling在Bayeisan中十分重要。比较高级的方法之一就是基于随机微分方程的抽样，例如这里的Langevin和HMC(Hamiltonian Monte Carlo)。

Langevin dynamic

下面是一个目标分布是标准正态分布的Langevin算法的小实验。分别是运行长度为100000、长度为1000000、更改噪声尺度得到的结果：

langevin的理论中，步长满足一定条件下，离散化导致的误差被加入的随机噪声掩盖，后期MH拒绝可能性yanhua会很低，几乎可忽略。
此类算法要求的计算量很不小。在分布容易模拟的情况下千万不要用。
增加噪声的尺度会直接影响最后的分布更加不集中。课件噪声在这个过程中起到很重要的作用。

Hamiltonian Monte Carlo

HMC领域一个重要的文章是Tianqi Chen大佬写的那篇Stochastic Gradient Hamiltonian Monte Carlo。

相比于引入摩擦这个概念，我觉得作者在离散化和SDE间的转换更有意思。固定了步长ϵ\epsilonϵ后，导致的误差本来是O(ϵ2)O(\epsilon^2)O(ϵ2)，直接连续化到SDE对应的并不是一个布朗运动噪声（对应的东西应该并不存在），所以需要更换一种方法。作者将O(ϵ2)O(\epsilon^2)O(ϵ2)的离散方差解释为布朗运动N(0,ϵVdt)N(0, \epsilon Vdt)N(0,ϵVdt)的ϵ\epsilonϵ离散。最终和Langevin结合起来，因此可以扔掉MH step.

上述内容有一个有趣的问题：连续化的SDE里的噪音一直带有ϵVdt\epsilon VdtϵVdt的方差，ϵ\epsilonϵ是前面的离散步长。

作者在不考虑ϵ\epsilonϵ为离散步长的前提下，继续对SDE做分析，加入同样带有ϵ\epsilonϵ的摩擦力项。发现得到的SDE是满足目标分布平稳的，再对之做ϵt\epsilon_tϵt的离散。

逻辑有点儿复杂，之后再想想合理性怎么解释。

# target: Normal(0,1)library(ggplot2)# Langevin
L = 100000
e = 0.001x <- c(0)
for(i in 1:L)
{x[i+1] = x[i] + e/2*(-x[i]) + sqrt(e)*rnorm(1)
}X = x[as.integer(L/2):L]
ggplot()+geom_histogram(aes(x=X,colour="Langevin dynamic",..density..),binwidth = 0.2)+geom_line(aes(x=seq(-3,3,0.1),y=dnorm(seq(-3,3,0.1)),colour="True density"))# more stepsL = 1000000
e = 0.001x <- c(0)
for(i in 1:L)
{x[i+1] = x[i] + e/2*(-x[i]) + sqrt(e)*rnorm(1)
}X = x[as.integer(L/2):L]
ggplot()+geom_histogram(aes(x=X,colour="Langevin dynamic",..density..),binwidth = 0.2)+geom_line(aes(x=seq(-3,3,0.1),y=dnorm(seq(-3,3,0.1)),colour="True density"))# double noise. Wrong result. L = 1000000
e = 0.001x <- c(0)
for(i in 1:L)
{x[i+1] = x[i] + e/2*(-x[i]) + 2*sqrt(e)*rnorm(1)
}X = x[as.integer(L/2):L]
ggplot()+geom_histogram(aes(x=X,colour="Wrong Langevin dynamic",..density..),binwidth = 0.2)+geom_line(aes(x=seq(-3,3,0.1),y=dnorm(seq(-3,3,0.1)),colour="True density"))