使用Markov Chain采样器表示采样分布- MATLAB & Simulink - MathWorks Deutschland

使用Markov Chain采样器表示采样分布

For more complex probability distributions, you might need more advanced methods for generating samples than the methods described in常见的伪随机数代生成方法。例如，在贝叶斯数据分析和大型组合问题中出现这种分布 Markov chain Monte Carlo (MCMC) simulations. An alternative is to construct a Markov chain with a stationary distribution equal to the target sampling distribution, using the states of the chain to generate random numbers after an initial burn-in period in which the state distribution converges to the target.

使用Metropolis-Hastings算法

Metropolis-Hastings算法从仅以常数已知的分布绘制样本。从具有概率密度函数的分布生成随机数，该概率密度函数等于或与提议功能成比例。

To generate random numbers:

假设初始值x(t）。
Draw a sample,y(t), from a proposal distributionq(y|x(t)).
接受y(t) as the next samplex(t+ 1）具有概率r(x(t),y(t）），并保持x(t) as the next samplex(t+ 1）具有概率1 -r(x(t),y(t)), where:

$r (x, y) = m i n {\frac{f (y)}{f (x)} \frac{q (x | y)}{q (y | x)}, 1}$
Incrementt→t+ 1, and repeat steps 2 and 3 until you get the desired number of samples.

Generate random numbers using the Metropolis-Hastings method with themhsamplefunction. To produce quality samples efficiently with the Metropolis-Hastings algorithm, it is crucial to select a good proposal distribution. If it is difficult to find an efficient proposal distribution, use slice sampling (片状）或汉密尔顿蒙特卡罗（HMCSampler.) instead.

使用切片抽样

In instances where it is difficult to find an efficient Metropolis-Hastings proposal distribution, the slice sampling algorithm does not require an explicit specification. The slice sampling algorithm draws samples from the region under the density function using a sequence of vertical and horizontal steps. First, it selects a height at random from 0 to the density functionf(x）。然后，它选择一个新的x通过从所选高度上方的密度的水平“切片”采样来随机测量值。类似的切片采样算法用于多变量分布。

If a functionf(x) proportional to the density function is given, then do the following to generate random numbers:

假设初始值x(t) within the domain off(x）。
Draw a real valueyuniformly from (0,f(x(t））），从而定义水平“切片”为S= {x:y<f(x)}.
找到一个间隔I= (L,R) aroundx(t) that contains all, or much of the “slice”S。
画出新的点x(t+ 1) within this interval.
Incrementt→t+ 1 and repeat steps 2 through 4 until you get the desired number of samples.

切片采样可以从具有任意形式的密度函数的分布生成随机数，只要有效的数值可以找到间隔I= (L,R), which is the “slice” of the density.

Generate random numbers using the slice sampling method with the片状function.

使用汉密尔顿蒙特卡洛

Metropolis-Hastings and slice sampling can produce MCMC chains that mix slowly and take a long time to converge to the stationary distribution, especially in medium-dimensional and high-dimensional problems. Use the gradient-based Hamiltonian Monte Carlo (HMC) sampler to speed up sampling in these situations.

To use HMC sampling, you must specify logf(x)（最多一种添加剂常数）及其梯度。您可以使用数值梯度，但这导致采样较慢。所有采样变量必须是不受约束的，这意味着日志f(x)它的渐变是全部真实定义的x。To sample constrained variables, transform these variables into unconstrained ones before using the HMC sampler.

The HMC sampling algorithm introduces a random “momentum vector”z并定义一个关节密度zand the “position vector”xasp（x，z）= f（x）g（z）。目标是从这个联合分布中进行采样，然后忽略这些值z— the marginal distribution ofx具有所需的密度f(x)。

HMC算法将高斯密度与协方差矩阵分配M（“质量矩阵”）到z:

$g (z) \propto \exp (- \frac{1}{2} z^{T} M^{- 1} z)$

Then, it defines an “energy function” as

$E (x, z) = - \log f (x) + \frac{1}{2} z^{T} M^{- 1} z = U (x) + K (z)$

withU(x) =- 日志f(x)“潜在能量”和k（z）= z^TM^-1z / 2.the “kinetic energy”. The joint density is given byP(x,z) ∝exp {-E(x,z)}。

要生成随机样本，HMC算法：

假设初始值xof the position vector.
生成动量矢量样本：z ∼ g(z)。
Evolves the state(x, z)for some amount of fictitious timeτ.to a new state(x’,z’)using the “equations of motion”:

$\frac{d z}{d τ.} = - \frac{\partial U}{\partial x}$

$\frac{d x}{d τ.} = \frac{\partial K}{\partial z}$

如果可以确切地解决运动方程，则能量（并且因此密度）将保持恒定：E（x，z）= e（x'，z'）。In practice, the equations of motions must be solved numerically (usually using so-called leapfrog integration) and the energy is not conserved.
接受x’as the next sample with probabilityp_acc= min(1, exp{E(x,z) – E(x’,z’)}），并保持x作为下一个具有概率1的样本 -p_acc。
重复步骤2到4，直到它产生所需的样本数量。

To use HMC sampling, create a sampler using theHMCSampler.function. After creating a sampler, you can compute MAP (maximum-a-posteriori) point estimates, tune the sampler, draw samples, and check convergence diagnostics. For an example of this workflow, see贝叶斯Linear Regression Using Hamiltonian Monte Carlo。