在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中截斷多元正態？

回答

相關問題