python - SciPy中的多元正态被截断了吗？

我正在尝试自动化某个过程，该过程有时需要从截断的多元法线中提取样本。也就是说，它是一个正态多元正态分布(即高斯分布)，但变量被约束为一个长方体。我给定的输入是完整多元正态的均值和协方差，但是我需要在框中输入样本。

到目前为止，我只是拒绝在盒子外面采样，并根据需要重新采样，但是我开始发现我的过程有时会给我(a)大协方差，并且(b)接近边缘。这两个事件共同影响了我的系统速度。

所以我想做的是首先正确地对分布进行采样。谷歌搜索仅导致this discussion或truncnorm中的 scipy.stats distribution。前者没有定论，而后者似乎是针对一个变量的。是否有任何本地多元截断法线？它会比拒绝 sample 更好，还是我应该做些更聪明的事情？

我将开始研究自己的解决方案，即将未截断的高斯旋转到主轴(使用SVD分解等)，使用截断的高斯积对分布进行采样，然后将其旋转回去，并根据需要拒绝/重新采样。如果截短的采样效率更高，我认为这应该更快地采样所需的分布。

最佳答案

因此，根据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通过将参数threads=Nthreads添加到EnsembleSampler初始化中轻松支持基本并行化。因此，您可以快速进行此操作。

关于python - SciPy中的多元正态被截断了吗？，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/20115917/