使用正態分佈的繪圖模擬均勻分佈的繪圖

October 2, 2014

我最近購買了一個數據科學面試資源，其中一個概率問題如下：

給定具有已知參數的正態分佈的繪製，如何模擬均勻分佈的繪製？

我最初的想法是，對於一個離散的隨機變量，我們可以將正態分佈分解為 K 個獨特的子部分，其中每個子部分在正態曲線下的面積相等。然後我們可以通過識別變量最終落入正態曲線的哪個區域來確定變量採用哪個 K 值。

但這僅適用於離散隨機變量。我對如何對連續隨機變量做同樣的事情進行了一些研究，但不幸的是，我只能找到像逆變換採樣這樣的技術，它可以將一個均勻隨機變量用作輸入，並且可以從其他分佈中*輸出隨機變量。*我在想也許我們可以反過來做這個過程來獲得統一的隨機變量？

我還考慮過可能使用 Normal 隨機變量作為線性同餘生成器的輸入，但我不確定這是否可行。

關於我如何處理這個問題的任何想法？

本著使用與計算正態分佈無關的簡單代數計算的精神，我傾向於以下內容。它們是按照我的想法排序的（因此需要變得越來越有創意），但我把最好的——也是最令人驚訝的——留到了最後。

反轉Box-Mueller 技術：從每對法線，兩個獨立的製服可以構造為（在區間）和（在區間）。

取兩個一組的法線並將它們的平方相加以獲得一個序列變量. 從對獲得的表達式

會有一個分佈，是均勻的。

這只需要基本的、簡單的算術應該是清楚的。 3. 因為來自標準雙變量正態分佈的四對樣本的 Pearson 相關係數的精確分佈均勻分佈在，我們可以簡單地取四對一組的法線（即每組有八個值）並返回這些對的相關係數。（這涉及簡單的算術加上兩個平方根運算。） 4. 自古以來就知道球體（三空間中的一個表面）的圓柱投影是等積的。這意味著在球體上均勻分佈的投影中，水平坐標（對應經度）和垂直坐標（對應緯度）都將具有均勻分佈。因為三元標準正態分佈是球對稱的，所以它在球面上的投影是均勻的。獲得經度本質上與 Box-Mueller 方法 ( qv ) 中的角度計算相同，但投影緯度是新的。到球體上的投影僅歸一化三坐標那時是投影緯度。因此，將 Normal 變量分成三組，, 併計算

為了. 5. 因為大多數計算系統以二進製表示數字，所以統一數字的生成通常從在它們之間生成均勻分佈的*整數開始。*和（或一些高功率與計算機字長有關）並根據需要重新調整它們。這樣的整數在內部表示為字符串二進制數字。我們可以通過將 Normal 變量與其中值進行比較來獲得獨立的隨機位。因此，只需將 Normal 變量分成大小等於所需位數的組，將每個變量與其平均值進行比較，並將結果的真/假結果序列組合成二進制數即可。寫作對於位數和對於符號（即，什麼時候和否則）我們可以將得到的歸一化統一值表示為用公式

變量可以從任何中位數為的連續分佈中得出（如標準法線）；它們被處理成組每組產生一個這樣的偽均勻值。 6. 拒絕抽樣是一種從任意分佈中抽取隨機變量的標準、靈活、強大的方法。假設目標分佈有 PDF. 一個值是按照另一種分佈用PDF繪製的. 在拒絕步驟中，一個統一的值介於兩者之間和獨立於和相比: 如果它更小，保留，否則重複該過程。不過，這種方法似乎是循環的：我們如何通過需要一個統一變量開始的過程生成一個統一變量？

答案是我們實際上不需要統一變量來執行拒絕步驟。相反（假設) 我們可以擲一枚公平的硬幣來獲得或者隨機。這將被解釋為統一變量的二進製表示中的第一位在區間. 當結果是, 這意味著; 否則，. *一半的時間，這足以決定拒絕步驟：*如果但硬幣是,應該被接受；如果但硬幣是,應該被拒絕；否則，我們需要再次擲硬幣以獲得下一位. 因為——無論什麼價值有——有一個每次翻轉後停止的機會，預期的翻轉次數僅為.

如果預期的拒絕數量很少，拒絕抽樣可能是值得的（且有效的）。我們可以通過在 Normal PDF 下擬合最大可能的矩形（表示均勻分佈）來實現這一點。

使用微積分優化矩形的面積，你會發現它的端點應該位於, 它的高度等於，使其面積略大於. 通過使用這個標準的法向密度作為並拒絕區間外的所有值自動地，否則應用拒絕程序，我們將獲得統一的變量有效率的：

在一小部分當時，Normal 變量超出了並立即被拒絕。(是標準的 Normal CDF。）

在剩下的部分時間裡，必須遵循二元拒絕程序，平均需要兩個以上的 Normal 變量。

整個過程需要平均腳步。產生每個統一結果所需的正態變量的預期數量計算為

雖然這非常有效，但請注意（1）計算 Normal PDF 需要計算指數和（2）值必須一勞永逸地預先計算。它的計算量仍然比 Box-Mueller 方法 ( qv ) 少一些。 7. 均勻分佈的順序統計有指數級的差距。由於兩個法線（均值為零）的平方和是指數的，我們可以生成一個實現獨立制服，通過對此類法線對的平方求和，計算這些法線的累積和，重新調整結果以落在區間內，並刪除最後一個（始終等於）。這是一種令人愉悅的方法，因為它只需要平方、求和和（最後）一個除法。

這值將自動按升序排列。如果需要這樣的排序，這種方法在計算上優於所有其他方法，因為它避免了某種成本。但是，如果需要一系列獨立的製服，則對這些制服進行排序隨機值可以解決問題。由於（如 Box-Mueller 方法中所見，qv）每對法線的比率獨立於每對的平方和，我們已經有了獲得隨機排列的方法：按相應的比率對累積和進行排序. （如果非常大，這個過程可以在較小的群體中進行幾乎沒有效率損失，因為每個組只需要要創建的法線統一的價值觀。對於固定，因此漸近計算成本為=, 需要要生成的正態變量統一的值。） 8. 一個極好的近似值，任何具有大標準偏差的正態變量在小得多的值的範圍內看起來都是一致的。將此分佈滾動到範圍內（通過只取值的小數部分），我們因此獲得了一個對於所有實際目的來說都是均勻的分佈。這是非常有效的，需要最簡單的算術運算之一：只需將每個 Normal 變量向下舍入到最接近的整數並保留多餘的部分。當我們檢查實際實現時**，這種方法的簡單性變得引人注目：**R
rnorm(n, sd=10) %% 1
可靠地產生n範圍內的統一值n僅以Normal 變量和幾乎沒有計算為代價。

（即使標準差為, 這個近似的 PDF 與一個統一的 PDF 不同，如下圖所示，相差不到一個！為了可靠地檢測它，需要一個樣本值——這已經超出了任何標準隨機性測試的能力。使用較大的標準偏差，不均勻性非常小，甚至無法計算。例如，SD 為如代碼所示，與統一 PDF 的最大偏差僅為.)

在每種情況下，“具有已知參數”的法線變量都可以很容易地重新定位並重新調整為上面假設的標準法線。之後，生成的均勻分佈的值可以重新定位並重新縮放以覆蓋任何所需的間隔。這些只需要基本的算術運算。

下面的R代碼證明了這些構造的簡單性，其中大多數只使用一兩行代碼。它們的正確性由基於每種情況下的獨立值（所有七個模擬都需要大約 12 秒）。作為參考 - 如果您擔心這些圖中出現的變化量 -R最後包含使用的統一隨機數生成器模擬的統一值的直方圖。

所有這些模擬均使用測試基於垃圾箱；沒有一個可以被認為是顯著不均勻的（最低的 p 值是–對於由R的實際統一數字生成器生成的結果！）。
set.seed(17)
n <- 1e5
y <- matrix(rnorm(floor(n/2)*2), nrow=2)
x <- c(atan2(y[2,], y[1,])/(2*pi) + 1/2, exp(-(y[1,]^2+y[2,]^2)/2))
hist(x, main="Box-Mueller")

y <- apply(array(rnorm(4*n), c(2,2,n)), c(3,2), function(z) sum(z^2))
x <- y[,2] / (y[,1]+y[,2])
hist(x, main="Beta")

x <- apply(array(rnorm(8*n), c(4,2,n)), 3, function(y) cor(y[,1], y[,2]))
hist(x, main="Correlation")

n.bits <- 32; x <-  (2^-(1:n.bits)) %*% matrix(rnorm(n*n.bits) > 0, n.bits)
hist(x, main="Binary")

y <- matrix(rnorm(n*3), 3)
x <- y[1, ] / sqrt(apply(y, 2, function(x) sum(x^2)))
hist(x, main="Equal area")

accept <- function(p) { # Using random normals, return TRUE with chance `p`
 p.bit <- x <- 0
 while(p.bit == x) {
   p.bit <- p >= 1/2
   x <- rnorm(1) >= 0
   p <- (2*p) %% 1
 }
 return(x == 0)
}
y <- rnorm(ceiling(n * sqrt(exp(1)*pi/2))) # This aims to produce `n` uniforms
y <- y[abs(y) < 1]
x <- y[sapply(y, function(x) accept(exp((x^2-1)/2)))]
hist(x, main="Rejection")

y <- matrix(rnorm(2*(n+1))^2, 2)
x <- cumsum(y)[seq(2, 2*(n+1), 2)]
x <- x[-(n+1)] / x[n+1]
x <- x[order(y[2,-(n+1)]/y[1,-(n+1)])] 
hist(x, main="Ordered")

x <- rnorm(n) %% 1 # (Use SD of 5 or greater in practice)
hist(x, main="Modular")

x <- runif(n)      # Reference distribution
hist(x, main="Uniform")