4

我想模拟来自非标准密度函数的数据。我已经找到了以下链接（如何最好地使用其概率函数模拟任意单变量随机变量？）。但是，这会产生奇怪的结果。不知何故，这个累积密度函数 ( 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 中的数据？

4

3 回答 3

5

正如@pjs 所指出的，我们可以使用拒绝采样（有关详细信息，请查看 wiki）。

这是这种方法的一种实现。

最重要的一步是找到一个分布 g，我们可以从中进行采样，并且它存在于 M 中，使得所有点的 M * g > f

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")

在此处输入图像描述

现在，我们可以执行算法了

set.seed(42)
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] 19.7086

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

在此处输入图像描述

录取率不是很高。您可以尝试找到更好的包络函数或使用 Metropolis Hasting 算法。

于 2013-04-21T21:04:32.210 回答
4

如果积分区间很大，那么密度的峰值很难找到：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-21T19:10:47.277 回答
0

当我玩弄你的 CDF 时，很快发现大部分动作都是针对 180 到 350 之间的 x 进行的，我通过绘制该范围内的密度来确认这一点。

我很确定 x = 9702 处的结果反映了当您涉及 25 次方和 26 次方时计算的数值不稳定性。如果您不信任您的 CDF 或者它不可逆，则另一个基于 pdf 的选项是接受/拒绝。您应该能够使用一个最小值 = 180、最大值约为 300、众数约为 200 的简单三角形作为边界函数 g(x)，并遵循 Wikipedia 上描述的算法以获得相当好的结果。

一般来说，如果反转不适用于任意分布，您的其他选择是 1）基于 pdf 相对于边界函数的接受/拒绝，2）组合（您可以将分布解构为更易于生成的组件并选择一个使用条件概率的适当组件），或 3）“特殊技巧” - 是否存在卷积或参数化给出分布等价的情况（例如，N(0,1)^2 = chi-square(1), chi-square(k) = k 个独立卡方 (1) 的总和，exp(2) = 卡方 (2) 等...)。请参阅Luc Devroye 的关于非均匀随机变量生成的书，以全面了解您的选择。

于 2013-04-21T19:41:57.660 回答