Генерация случайных чисел Ломакса в R - PullRequest
Новая
       33

Генерация случайных чисел Ломакса в R

0 голосов
/ 03 июля 2018

Как я могу генерировать случайные числа Lomax (тип Паретто II), используя R?

Если U∈ [0,1) является равномерно распределенной случайной величиной, то

л (хт, α) = Р (хт, α) -xm

генерирует распределенную случайную величину Lomax.

1 Ответ

0 голосов
/ 04 июля 2018

В качестве альтернативы использованию VGAM::rlomax нетрудно написать свой выигранный Lomax генератор случайных чисел с использованием выборки обратного преобразования .

Cdf распределения Lomax с параметрами shape и scale задается как F(x) = 1 - (1 + x / scale)^(-alpha). Все, что нам нужно сделать, это решить F(F^(-1)(x)) = x для F^(-1)(x), где x ~ Unif(0, 1).

С этим решением мы можем определить следующую функцию для рисования случайных выборок Lomax 

rlomax.its <- function(N, scale, shape) {
    scale * ((1 - runif(N)) ^ (-1/shape) - 1)
}

Теперь мы берем N = 1e5 выборок из распределения Lomax с scale = 1 и shape = 2 и сравниваем с выборками из VGAM::rlomax

library(VGAM);
N <- 1e5;
set.seed(2017);
x.VGAM <- rlomax(N, scale = 1, shape3.q = 2)
x.ITS <- rlomax.its(N, scale = 1, shape = 2)

summary(x.VGAM);
#Min.  1st Qu.   Median     Mean  3rd Qu.     Max.
#0.0000   0.1536   0.4143   0.9985   1.0006 925.0784

summary(x.ITS);
#Min.  1st Qu.   Median     Mean  3rd Qu.     Max.
#0.0000   0.1548   0.4158   1.0016   1.0086 280.3248

Давайте сравним графики плотности для разных размеров выборки, используя оба метода. 

set.seed(2017);
bind_rows(map(
    setNames(2:5, paste0("N=10^", 2:5)),
    ~list(ITS = rlomax.its(10^(.x), 1, 2), VGAM = rlomax(10^(.x), 1, 2))),
    .id = "N") %>%
    gather(key, value, -N) %>%
    ggplot(aes(log10(value), fill = key)) +
    geom_density(alpha = 0.4) +
    facet_wrap(~ N)

enter image description here

Очевидно, что по мере увеличения * 1034 распределения обоих методов сходятся.

Относительно того, какой метод быстрее, мы можем запустить быстрое microbenchmark на основе рисования N=1e6 семплов Lomax от обоих методов 

library(microbenchmark);
res <- microbenchmark(
    ITS = rlomax.its(1e6, 1, 2),
    VGAM = rlomax(1e6, 1, 2))
#Unit: milliseconds
# expr       min        lq      mean    median        uq      max neval cld
#  ITS  79.22709  84.11703  88.48358  86.29181  91.07074 109.3536   100  a
# VGAM 159.56578 175.88731 218.92212 183.09769 222.64697 359.9311   100   b

library(tidyverse)
autoplot(res)

enter image description here

Давайте посмотрим на зависимость времени выполнения как функции отрисованных сэмплов 

library(tidyverse);
library(ggthemes);
res <- map_df(seq(2, 6, length.out = 20), function(x)
    cbind(x = 10^(x), microbenchmark(
        ITS = rlomax.its(10^(x), 1, 2),
        VGAM = rlomax(10^(x), 1, 2))))
res %>%
    mutate(N = factor(as.numeric(factor(x)))) %>%
    ggplot(aes(x = N, y = log10(time), colour = expr)) +
    geom_tufteboxplot(outlier.colour="transparent") +
    theme_minimal() +
    scale_x_discrete(
        breaks = c(1, 5, 10, 15, 20),
        labels = paste0("10^", 2:6))

enter image description here

У меня нет времени, чтобы исследовать это дальше, но оказывается, что в среднем метод обратной выборки немного (но последовательно) быстрее.

...