Probability

在單位圓和單位正方形之間高效生成點

  • January 26, 2017

我想從此處定義的藍色區域生成樣本:

在此處輸入圖像描述

天真的解決方案是在單位平方中使用拒絕抽樣,但這僅提供了一個(~21.4%) 效率。

有什麼方法可以更有效地採樣嗎?

每秒200萬點會做嗎?

分佈是對稱的:我們只需要計算出整個圓的八分之一的分佈,然後將其複製到其他八分圓。在極坐標中, 角度的累積分佈對於隨機位置在價值由三角形之間的面積給出和從延伸的圓弧到. 因此它與

它的密度從何而來

我們可以從這個密度中取樣,比如說,一種拒絕方法(它有效率).

徑向坐標的條件密度正比於之間和. 這可以通過 CDF 的簡單反演進行採樣。

如果我們生成獨立樣本, 轉換回笛卡爾坐標採樣這個八分圓。因為樣本是獨立的,所以隨機交換坐標會根據需要從第一象限產生一個獨立的隨機樣本。(隨機交換只需要生成一個二項式變量來確定要交換多少個實現。)

每一個這樣的實現平均而言,需要一個統一變量(對於) 更多的乘以兩個均勻變量(對於) 和少量(快速)計算。那是每個點變化(當然,有兩個坐標)。完整的詳細信息在下面的代碼示例中。這個數字繪製了超過 50 萬個點中的 10,000 個。

數字

這是R生成此模擬並對其進行計時的代碼。

n.sim <- 1e6
x.time <- system.time({
 # Generate trial angles `theta`
 theta <- sqrt(runif(n.sim)) * pi/4
 # Rejection step.
 theta <- theta[runif(n.sim) * 4 * theta <= pi * tan(theta)^2]
 # Generate radial coordinates `r`.
 n <- length(theta)
 r <- sqrt(1 + runif(n) * tan(theta)^2)
 # Convert to Cartesian coordinates.
 # (The products will generate a full circle)
 x <- r * cos(theta) #* c(1,1,-1,-1)
 y <- r * sin(theta) #* c(1,-1,1,-1)
 # Swap approximately half the coordinates.
 k <- rbinom(1, n, 1/2)
 if (k > 0) {
   z <- y[1:k]
   y[1:k] <- x[1:k]
   x[1:k] <- z
 }
})
message(signif(x.time[3] * 1e6/n, 2), " seconds per million points.")
#
# Plot the result to confirm.
#
plot(c(0,1), c(0,1), type="n", bty="n", asp=1, xlab="x", ylab="y")
rect(-1, -1, 1, 1, col="White", border="#00000040")
m <- sample.int(n, min(n, 1e4))
points(x[m],y[m], pch=19, cex=1/2, col="#0000e010")

引用自:https://stats.stackexchange.com/questions/258303

comments powered by Disqus