在SciPy中截断多元正态？

我试图自动化一个过程，在某些时候需要从截断的多元正态图中抽取样本。也就是说，它是一个正态多元正态分布（即高斯），但变量被限制为一个长方体。我给出的输入是完整多变量正态分布的均值和协方差，但我需要在我的框中输入样本。到目前为止，我只是拒绝了盒子外的样品并根据需要重新采样，但我开始发现我的过程有时会给我（a）大的协方差和（b）意味着接近于边缘。这两个事件与我的系统速度相悖。在SciPy中截断多元正态？

所以我想要做的是首先正确地分配分配。谷歌搜索仅导致this discussion或truncnorm distribution在scipy.stats。前者不确定，后者似乎是一个变量。是否有任何原生多变量截断正态？它会比拒绝样品更好，还是我应该做更聪明的事情？

我打算开始研究自己的解决方案，它将未经截断的高斯旋转到它的主轴（用SVD分解或某物），使用截断高斯的乘积来对分布进行采样，然后旋转取样返回，并根据需要拒绝/重新取样。如果截断采样更有效，我认为这应该更快地对期望的分布进行采样。

来源

2013-11-21 Warrick

因此，根据the Wikipedia article，采样多元截断正态分布（MTND）更困难。我最终采取了一种相对简单的方法，并使用MCMC采样器来放松对MTND的初始猜测，如下所示。

我用emcee来做MCMC的工作。我发现这个软件包非常易于使用。它只需要一个返回所需分布的对数概率的函数。所以，我这个函数

from numpy.linalg import inv 

def lnprob_trunc_norm(x, mean, bounds, C): 
    if np.any(x < bounds[:,0]) or np.any(x > bounds[:,1]): 
     return -np.inf 
    else: 
     return -0.5*(x-mean).dot(inv(C)).dot(x-mean)

这里，C是多元正态分布的协方差矩阵。然后，您可以运行类似

S = emcee.EnsembleSampler(Nwalkers, Ndim, lnprob_trunc_norm, args = (mean, bounds, C)) 

pos, prob, state = S.run_mcmc(pos, Nsteps)

给定mean，bounds和C。您需要为步行者的立场pos初始猜测，这可能是围绕平均值球，

pos = emcee.utils.sample_ball(mean, np.sqrt(np.diag(C)), size=Nwalkers)

或由未截断多元正态分布采样，

pos = numpy.random.multivariate_normal(mean, C, size=Nwalkers)

等。我个人首先做了几千个样本丢弃步骤，因为它很快，然后强制剩余的异常值回到范围内，然后运行MCMC采样。

收敛的步骤数由您决定。

还要注意，emcee通过在EnsembleSampler初始化中添加参数threads=Nthreads可轻松支持基本并行化。所以你可以快速地让这个炽热的。

来源

2013-11-25 09:14:33 Warrick

在SciPy中截断多元正态？

回答

相关问题