Как я мог воспроизвести этот сюжет Log-log для распределения Pareto в r?

Question

Я пытаюсь повторить этот сюжет

, который, согласно статье, построил полную линию, генерируя случайные числа с этимуравнение

Код, который я использовал для генерации случайных чисел:

r <- c(runif(600, min = 0, max = 1))
pnumbers = c()
a = 0.17
b = 4200
for(i in 1:600){
  pnumbers[i] = a*(1 - r[i]*(1 - (a/b)^0.75))^(-1/0.75)
}
x2 <- sort(pnumbers)

и использует эти случайные числа в уравнении11 показано на этом рисунке

Эта функция была рассчитана с помощью этого кода

pareto1 <- ppareto(x2, 0.17, 0.75)
survpareto <- 1 - pareto1

Я мог получить прямую пунктирную линию, ноЯ не смог получить график кривой с пустыми кругами.Предполагается, что график пустых циклов был сделан из набора данных и уравнения 11, но я получил еще одну прямую линию!На самом деле ... та же самая прямая линия! График, который я получил Тот график, который я получил с этим кодом

pdf("PruebaGraficoLogLog.pdf")
pareto2 <- 1 - (0.17/x)^0.75
survpareto2 <- 1 - pareto2
plot(x2, survpareto, log = "xy", col = "blue", type = "l", lty = 5)
points(x, survpareto2, log = "xy")
dev.off()

Мой вопрос: что мне нужно сделать, чтобы правильно воспроизвести первый план?Что я делаю неправильно?Спасибо за вашу помощь и сотрудничество.

РЕДАКТИРОВАТЬ: Я изменил название, чтобы сделать его более конкретным и подробным.Это преамбула (пакеты), которую я использовал в своей игрушечной модели

library(EnvStats)
library(stats)
library(base)

Я не использовал пакет fitdistrplus, потому что (и я не знаю почему) я не смог установить его вмой компьютер.Я использую R 3.3, но установка пакета всегда заканчивается неудачей.

Answer 1

Сначала два важных комментария:

• Распределение выживаемости ваших выборок, показанное кружками на рисунке, делает не параметры соответствия a, b, c, которые вы даете в своем посте.Можете ли вы объяснить, как вы пришли к этим конкретным значениям: a = 0.17, b = 4200, c = 0.75?

• Я должен сказать, что я не совсем понимаю весь смыслупражнения.На рисунке показаны выборки из ограниченного / усеченного распределения Парето (построены с использованием выборки с обратным преобразованием, см. Ниже);затем показано, что распределение выживаемости выборок согласуется с распределением ограниченного распределения Парето, а не неограниченным распределением Парето (что, очевидно, соответствует ожиданиям).Обычно это то, что вы делаете, когда пытаетесь оценить параметров из базового распределения (здесь: ограниченный Парето). Так, возможно, вы спрашиваете, как оценить параметры усеченного распределения Парето? Если это так, это будет зависеть от данных ( в случае случайных данных, нет фиксированного начального числаозначает отсутствие воспроизводимости ) и метод оценки (обычно ML).Возможно, будет полезно взглянуть на MASS:fitdistr.

---

Помимо этих комментариев, здесь приведен воспроизводимый пример для генерации x и построения log x против log S(x).

1. Создайте n=600 образцов x в соответствии с уравнением, которое вы даете (я предполагаю, что это уравнение 13 из подписи к рисунку).

   set.seed(2018)
   rsample <- function(n, a, b, c) a * (1 - runif(n) * (1 - (a / b)^c)) ^ (-1 / c)
   
   x <- rsample(600, 0.17, 4200, 0.75)
   

   Обратите внимание, что x генерируется с помощью выборки с обратным преобразованием (ITS) из ограниченного распределения Парето .Легко отобразить коэффициенты a, b, c на коэффициенты из определения ограниченного распределения Парето из Википедии:

   a = L         (location parameter)
   b = H         (location parameter)
   c = alpha     (shape parameter)
   

2. Мы вычисляем эмпирическое кумулятивное распределениефункция F_X(x) = P(X ≤ x) с использованием ecdf

   Px <- ecdf(x)
   

3. Теперь мы можем вычислить F_X(x) для любых значений x (учитывая, что поддержка ограниченного распределения Парето равна L ≤ x ≤ H).Мы выбираем значения так, чтобы они соответствовали интервалу, показанному в лог-масштабе на рисунке.Тогда функция выживания просто S_X(x) = 1 - F_X(x) = P(X > x).

   library(tidyverse)
   df <- data.frame(
       x = exp(seq(0, 1.6, 0.05))) %>%
       mutate(
           Px = Px(x),
           Sx = 1 - Px)
   

4. Мы строим x против функции выживания S_X(x) в логарифмическом масштабе.

   ggplot(df, aes(log(x), log(Sx))) +
       geom_point(size = 3, shape = 21)
   

Как видно, форма распределения выживания не соответствует форме на рисунке, демонстрируя, чтопараметры a = 0.17, b = 4200, c = 0.75 не согласуются с параметрами, используемыми для выборок функции выживания из рисунка.