Как генерировать случайные данные из функции плотности Парето в R?

Question

Я пытаюсь сгенерировать данные из заданной плотности Парето в R.

Плотность Парето: F (x) = | X | ^ (- 3) * 1 | x |> 1

Я знаю, что мне нужно использовать функцию rpareto из библиотеки актуаров, но я не уверен, как мне преобразовать заданную плотность парето в параметры.

Answer 1

Используя метод обратной выборки, формулы из Распределение Парето

Здесь обновленная версия с вычисленным средним и sd

rpar <- function(n, xm, a) {
    v <- runif(n)
    xm / v^(1.0/a)
}

rpar_mean <- function(xm, a) {
    result <- 1/0 # Inf
    if (a > 1.0)
        result <- a*xm/(a - 1.0)
    result
}

rpar_var <- function(xm, a) {
    result <- 1/0 # Inf
    if (a > 2.0)
        result <- xm*xm*a/((a - 1.0)^2*(a - 2.0))
    result
}

set.seed(54122345)

xm = 1.0
a  = 3.0

q <- rpar(10000, xm, a)
print(mean(q))
print(rpar_mean(xm, a))

print(sd(q))
print(sqrt(rpar_var(xm, a)))

Answer 2

Распределение Парето достаточно простое, чтобы учесть метод преобразования (он же Выборка обратного преобразования ).CDF распределения Парето:

$$ F (x) = \ int_ {x_ {min}} ^ xf (x ') \, dx' = 1- \ left (\ frac {x_ {min}} {x} \ right ^ k $$

, который можно аналитически преобразовать в

$$ F ^ {- 1} (y) = x_ {min} \ cdot (1-y) ^ {- 1 / k} $$

Таким образом, если вы преобразуете случайные числа, равномерно распределенные по $ (0,1) $, с помощью $ F ^ {- 1} $, вы получаете случайные числа, распределенные по Парето..

Редактировать: Извините, очевидно, латексный код здесь не поддерживается. Вот для вашего удобства код R:

k <- 5      # parameter k of the Pareto distribution
x.min <- 2  # cutoff point of Pareto distribution
N <- 500    # number of random points
x.random <- x.min*(1-runif(N))^(-1/k)

А вотпрактическая демонстрация того, что это работает:

h <- hist(x.random, freq=FALSE, plot=TRUE)
x <- seq(x.min, h$breaks[length(h$breaks)], by=0.01)
lines(x, k*x.min^k/x^(k+1), col="red")

график для k = 5 и x.min = 2