r - 如何在R中找到逆采样方法的逆

Question

通常对于逆采样方法，我们有一个密度，我们想从中采样。第一步是找到密度的累积密度函数。然后找到它的反函数，最后从均匀分布中找到随机采样值的反函数。

例如，我有这个函数y= ((3/2)/(1+x)^2)，所以 cdf 等于(3x)/2(x+1)并且 cdf 的倒数是((3/2)*u)/(1-(3/2)*u)

为了在 R 中做到这一点，我写了

 f<-function(x){
 y= ((3/2)/(1+x)^2)
 return(y)
}



cdf <- function(x){
  integrate(f, -Inf, x)$value
}

invcdf <- function(q){
  uniroot(function(x){cdf(x) - q}, range(x))$root
}
U <- runif(1e6)
X <- invcdf(U)

我有两个问题！首先：代码返回函数而不是样本。第二：有没有另一种简单的方法来完成这项工作？例如以更简单的方式找到 cdf 和逆？

我想补充一点，我不是在寻找代码的效率。我只是对初学者可以编写的代码感兴趣。

score 0 · Accepted Answer

您可以尝试使用数值方法进行逆采样。根据您的要求，这更多是关于方法的透明度而不是效率。

此函数将在给定范围内对给定函数进行数值积分（尽管它会修剪无限值）

cdf <- function(f, lower_bound, upper_bound)
{
  if(lower_bound < -10000) lower_bound <- -10000          # Trim large negatives
  if(upper_bound > 10000) upper_bound <- 10000            # Trim large positive
  x <- seq(lower_bound, upper_bound, length.out = 100001) # Finely divide x axis
  delta <- mean(diff(x))                                  # Get delta x (i.e. dx)
  mid_x <- (x[-1] + x[-length(x)])/2                      # Get the mid point of each slice
  result <- cumsum(delta * f(mid_x))                      # sum f(x) dx
  result <- result / max(result)                          # normalize
  list(x = mid_x, cdf = result)                           # return both x and f(x) in list
}

为了得到逆，我们在从 0 到 1 之间的均匀分布中抽取的随机数的 cdf 中找到最接近的值。然后我们看到 x 的哪个值对应于 cdf 的那个值。我们希望能够一次对 n 个样本执行此操作，因此我们使用sapply：

inverse_sample <- function(f, n = 1, lower_bound = -1000, upper_bound = 1000)
{
  CDF <- cdf(f, lower_bound, upper_bound)
  samples <- runif(n)
  sapply(samples, function(s) CDF$x[which.min(abs(s - CDF$cdf))])
}

我们可以通过绘制结果的直方图来测试它。我们将从正态分布的密度函数（dnorm在 R 中）开始，抽取 1000 个样本并绘制它们的分布：

hist(inv_sample(dnorm, 1000))

我们可以对指数分布做同样的事情，这次将积分限制设置在 0 到 100 之间：

hist(inv_sample(dexp, 1000, 0, 100))

最后我们可以用你自己的例子做同样的事情：

f <- function(x) 3/2/(1 + x)^2

hist(inv_sample(f, 1000, 0, 10))

r - 如何在R中找到逆采样方法的逆

1 回答 1

Related

Reference