Как найти обратное для метода обратной выборки в 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 и inverse более простыми способами?

Я хотел бы добавить, что я не ищу эффективность кода. Мне просто интересен код, который может быть написан новичком.

Answer 1

Вы можете попробовать численный подход к обратной выборке. Согласно вашему запросу, это больше касается прозрачности метода, чем эффективности.

Эта функция будет численно интегрировать данную функцию в заданный диапазон (хотя она будет обрезать бесконечные значения)

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
}

И чтобы получить обратное, мы находим ближайшее значение в cdf случайного числа, полученного из равномерного распределения между 0 и 1. Затем мы видим, какое значение 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))