模擬來自（非標準）密度函數的數據

我想模擬來自非標準密度函數的數據。我已經找到以下鏈接（How do I best simulate an arbitrary univariate random variate using its probability function?）。但是，這會產生奇怪的結果。不知何故，這個累積密度函數（cdf（））不能很好地工作。從一些價值觀來看，它給出了很奇怪的結果。例如，看看下面的代碼：模擬來自（非標準）密度函數的數據

density=function(x)(25*200.7341^25/x^26*exp(-(200.7341/x)^25)) 
cdf<-function(x) integrate(density,1,x)[[1]] 

cdf(9701) 
[1] 1 

cdf(9702) 
[1] 6.33897e-05

所以我的問題，我如何創建一個「好」的CDF功能？或者更直接地說，我如何模擬PDF中的數據？

來源

2013-04-21 Sjoerd Glaser

正如@pjs指出的那樣，我們可以使用Rejection sampling（請查看維基瞭解詳情）。

這是這種方法的一個實現。

最重要的一步是要找到一個分佈克從中我們可以採樣並從它存在M使得該M *克對所有點

f <- function(x) (25 * 200.7341^25/x^26 * exp(-(200.7341/x)^25)) 
g <- function(x) dnorm(x, mean = 200.7341, sd = 40) 
M <- 5 
curve(f, 0, 500) 
curve(M * g(x), 0, 500, add = TRUE, lty = "dashed")

enter image description here

>˚F現在，我們可以執行算法

k <- 1 
count <- 0 
res <- vector(mode = "numeric", length = 1000) 
while(k < 1001) { 
      z <- rnorm(n = 1, mean = 200.7341, sd = 40) 
      R <- f(z)/M * g(z) 
      if (R > runif(1)) { 
       res[k] <- z 
       k <- k + 1 
      } 
      count <- count + 1 
    } 

(accept_rate <- (k/count) * 100) 
## [1] 0.20773 

require(MASS) ## for truehist 
truehist(res) 
curve(f, 0, 250, add = TRUE)

enter image description here

驗收率很差，因此可以嘗試尋找更好的包絡函數或使用Metropolis Hasting算法。

來源

2013-04-21 21:04:32 dickoa

謝謝！這工作完美！雖然確實有點慢，但對於我想要做的事情來說，這很好。 – 2013-04-22 07:30:23

如果積分區間非常大，密度的峯值是很難找到：integrate可以輕易錯過它，，認爲要集成的功能是（幾乎）處處爲零。

如果您知道峯的位置，可以將峯積分前後的積分分爲三部分：。

# Density 
A <- 200.7341 
f <- function(x) 25*A^25/x^26 * exp(-(A/x)^25) 
a <- 150 
b <- 400 

# Numeric integration 
F1 <- function(x) { 
    if(x < a)  integrate(f, 1, x)[[1]] 
    else if(x < b) integrate(f, 1, a)[[1]] + integrate(f, a, x)[[1]] 
    else    integrate(f, 1, a)[[1]] + integrate(f, a, b)[[1]] + integrate(f, b, x)[[1]] 
} 

# Compare with the actual values 
F2 <- function(x) exp(-(A/x)^25) 
F1(200); F2(200) 
F1(1e4); F2(1e4) 
F1(1e5); F2(1e5) # Imprecise if b is too low...

確認您的間隔足夠大後，您可以刪除「之前」和「之後」間隔：其貢獻爲零。

F1 <- function(x) { 
    if(x < a)  0 
    else if(x < b) integrate(f, a, x)[[1]] 
    else    1 
}

來源

2013-04-21 19:10:47

當我與你的CDF發揮它周圍的迅速崛起，大部分的動作是180和350之間的x，我證實了繪製了該範圍內的密度。

我很確定x = 9702時的結果反映了當你有第25和第26次權力時計算的數值不穩定性。如果您不相信您的CDF或它不可逆，另一個基於pdf的選項是acceptance/rejection。你應該能夠使用一個簡單的三角形，最小值= 180，最大值300左右，模式200左右作爲邊界函數g（x），並遵循維基百科上描述的算法獲得相當好的結果。

一般而言，如果反演不適用於任意分佈，則其他選擇爲1）基於pdf的相對於邊界函數的接受/拒絕，2）合成（您是否可以將分佈解構爲更容易生成使用條件概率選擇一個合適的分量），或3）「特殊技巧」 - 是否存在卷積或參數化給出分佈等價的情況（例如，N（0,1）^ 2 =卡方（1）平方（k）= k獨立卡方（1）的總和，exp（2）=卡方（2）等）。請參閱關於非均勻隨機變量生成的Luc Devroye's book，以全面處理您的選項。

來源

2013-04-21 19:41:57 pjs

模擬來自（非標準）密度函數的數據

回答

相關問題