這個對嗎 ?(生成截斷範數多元高斯)
如果 IE,
我想要一個多變量情況下截斷正態分佈的類似版本。
更準確地說,我想生成一個規範約束(到一個值) 多元高斯英石
在哪裡
現在我觀察到以下內容:
如果,
因此通過選擇作為高斯樣本,可以限制作為截斷正態分佈的樣本(遵循高斯尾) 分配, 除了它的符號以概率隨機選擇.
現在我的問題是,
如果我生成每個向量樣本的作為,
和
在哪裡,,, (即截斷標量法線 RV
將要是一個範數約束 () 多元高斯?(即與定義如上)。我應該如何驗證?如果不是這樣,還有其他建議嗎?
編輯:
這是二維情況下點的散點圖,範數被截斷為高於“1”的值
注意:下面有一些很好的答案,但是沒有說明為什麼這個提議是錯誤的。事實上,這是這個問題的重點。
的多元正態分佈是球對稱的。您尋求的分佈會截斷半徑下面在. 因為這個標準只取決於長度,截斷分佈保持球對稱。自從與球面角無關和有一個distribution,因此您可以通過幾個簡單的步驟從截斷的分佈中生成值:
- 產生.
- 產生作為 a 的平方根分佈截斷在.
- 讓.
在第 1 步中,是作為一個序列獲得的標準正態變量的獨立實現。
在第 2 步中,很容易通過反轉分位數函數生成一個分佈:生成統一變量在(分位數)範圍內支持和並設置.
這是一個直方圖這種獨立的實現為了在尺寸,在下面截斷. 生成大約需要一秒鐘,證明了算法的效率。
紅色曲線是截斷的密度分佈按比例縮放. 它與直方圖的緊密匹配證明了這種技術的有效性。
要獲得截斷的直覺,請考慮以下情況,在方面。這是一個散點圖反對(為了獨立實現)。它清楚地顯示了半徑處的孔:
最後,注意(1)組件必須具有相同的分佈(由於球對稱)和(2),除非當,該共同分佈不是正態分佈。事實上,作為變大,(單變量)正態分佈的快速減小導致球面截斷的多元正態分佈的大部分概率聚集在表面附近-球體(半徑)。因此,邊際分佈必須接近一個縮放的對稱 Beta分佈集中在區間. 這在前面的散點圖中很明顯,其中在二維上已經很大了:點 limn a ring (a-sphere) 的半徑.
這是來自大小模擬的邊緣分佈的直方圖在尺寸與,(其中近似的 Beta分佈均勻):
自從第一次問題中描述的程序的邊緣是正常的(按構造),該程序不可能是正確的。
以下
R
代碼生成了第一個圖。它被構建為並行步驟 1-3 以生成.a
它被修改為通過更改變量,d
,來生成第二個圖形n
,sigma
然後在生成plot(y[1,], y[2,], pch=16, cex=1/2, col="#00000010")
後發出 plot 命令y
。一代在代碼中修改以獲得更高的數值分辨率:代碼實際生成並用它來計算.
根據假定的算法模擬數據、用直方圖對其進行總結和疊加直方圖的相同技術可用於測試問題中描述的方法。它將確認該方法無法按預期工作。
a <- 7 # Lower threshold d <- 11 # Dimensions n <- 1e5 # Sample size sigma <- 3 # Original SD # # The algorithm. # set.seed(17) u.max <- pchisq((a/sigma)^2, d, lower.tail=FALSE) if (u.max == 0) stop("The threshold is too large.") u <- runif(n, 0, u.max) rho <- sigma * sqrt(qchisq(u, d, lower.tail=FALSE)) x <- matrix(rnorm(n*d, 0, 1), ncol=d) y <- t(x * rho / apply(x, 1, function(y) sqrt(sum(y*y)))) # # Draw histograms of the marginal distributions. # h <- function(z) { s <- sd(z) hist(z, freq=FALSE, ylim=c(0, 1/sqrt(2*pi*s^2)), main="Marginal Histogram", sub="Best Normal Fit Superimposed") curve(dnorm(x, mean(z), s), add=TRUE, lwd=2, col="Red") } par(mfrow=c(1, min(d, 4))) invisible(apply(y, 1, h)) # # Draw a nice histogram of the distances. # #plot(y[1,], y[2,], pch=16, cex=1/2, col="#00000010") # For figure 2 rho.max <- min(qchisq(1 - 0.001*pchisq(a/sigma, d, lower.tail=FALSE), d)*sigma, max(rho), na.rm=TRUE) k <- ceiling(rho.max/a) hist(rho, freq=FALSE, xlim=c(0, rho.max), breaks=seq(0, max(rho)+a, by=a/ceiling(50/k))) # # Superimpose the theoretical distribution. # dchi <- function(x, d) { exp((d-1)*log(x) + (1-d/2)*log(2) - x^2/2 - lgamma(d/2)) } curve((x >= a)*dchi(x/sigma, d) / (1-pchisq((a/sigma)^2, d))/sigma, add=TRUE, lwd=2, col="Red", n=257)