Найти способ симулировать случайные числа для этого распределения

20

Я пытаюсь написать программу на R, которая имитирует псевдослучайные числа из распределения с помощью кумулятивной функции распределения:

F(Икс)знак равно1-ехр(-aИкс-бп+1Иксп+1),Икс0

гдеa,б>0,п(0,1)

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

Себастьян
источник
1
Не хватает времени для полного ответа, но вы можете проверить алгоритмы выборки по важности, как альтернативу.
Шеф
1
это не упражнение из учебника, я оговорил ограничение только потому, что это разумное предположение для моих данных
Себастьян
6
Затем меня удивляет «чудесная» нормализация с помощью (п+1)-1 которая превращает распределение в совершенную степень экспоненциального, но чудеса случаются (с небольшой вероятностью).
Сиань

Ответы:

49

Существует простое (и, если я могу добавить, элегантное) решение этого упражнения: поскольку 1F(x) выглядит как произведение двух распределений выживания:

(1-F(Икс))знак равноехр{-aИкс-бп+1Иксп+1}знак равноехр{-aИкс}1-F1(Икс)ехр{-бп+1Иксп+1}1-F2(Икс)
распределение F - это распределение
Иксзнак равномин{Икс1,Икс2}Икс1~F1,Икс2~F2
В этом случае F1 является экспоненциальным распределением Е(a) а F2 является 1/(п+1) -ым степень экспоненциального распределения Е(б/(п+1)) .

Соответствующий код R настолько прост, насколько это возможно

x=pmin(rexp(n,a),rexp(n,b/(p+1))^(1/(p+1))) #simulating an n-sample

и это определенно намного быстрее, чем обратное разрешение pdf и accept-reject:

> n=1e6
> system.time(results <- Vectorize(simulate,"prob")(runif(n)))
utilisateur     système      écoulé 
    89.060       0.072      89.124 
> system.time(x <- simuF(n,1,2,3))
utilisateur     système      écoulé 
     1.080       0.020       1.103 
> system.time(x <- pmin(rexp(n,a),rexp(n,b/(p+1))^(1/(p+1))))
utilisateur     système      écoulé 
     0.160       0.000       0.163 

с удивительно идеальной посадкой:

введите описание изображения здесь

Сиань
источник
5
действительно классное решение!
Себастьян
14

Вы всегда можете численно решить обратное преобразование.

Ниже я делаю очень простой поиск пополам. Для заданной входной вероятности (я использую так как в вашей формуле уже есть ), я начну с и . Затем я удваиваю до пор, пока . Наконец, я итеративно пополам интервал до тех пор, пока его длина не станет меньше а его средняя точка удовлетворит .QQпИксLзнак равно0Иксрзнак равно1ИксрF(Икср)>Q[ИксL,Икср]εИксMF(ИксM)Q

ECDF подходит вашему достаточно хорошо для моих вариантов и , и это достаточно быстро. Вероятно, вы могли бы ускорить это, используя некоторую оптимизацию типа Ньютона вместо простого поиска по разделению пополам.Faб

aa <- 2
bb <- 1
pp <- 0.1

cdf <- function(x) 1-exp(-aa*x-bb*x^(pp+1)/(pp+1))

simulate <- function(prob,epsilon=1e-5) {
    left <- 0
    right <- 1
    while ( cdf(right) < prob ) right <- 2*right

    while ( right-left>epsilon ) {
        middle <- mean(c(left,right))
        value_middle <- cdf(middle)
        if ( value_middle < prob ) left <- middle else right <- middle
    }

    mean(c(left,right))
}

set.seed(1)
results <- Vectorize(simulate,"prob")(runif(10000))
hist(results)

xx <- seq(0,max(results),by=.01)
plot(ecdf(results))
lines(xx,cdf(xx),col="red")

ECDF

С. Коласса - Восстановить Монику
источник
10

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

е(Икс)знак равно(a+бИксп)ехр{-aИкс-бп+1Иксп+1}
е(Икс)знак равноaе-aИксе-бИксп+1/(п+1)1+бИкспе-бИксп+1/(п+1)е-aИкс1
е(Икс)грамм(Икс)знак равноaе-aИкс+бИкспе-бИксп+1/(п+1)
граммξзнак равноИксп+1Иксзнак равноξ1/(п+1)
dИксdξзнак равно1п+1ξ1п+1-1знак равно1п+1ξ-пп+1
Иксимеет плотность вида где - нормализующая константа, тогда имеет плотность что означает , что (я) является распределяется как экспоненциальная переменная и (ii) константа равна единице. Следовательно, оказывается равным взвешенной смеси экспоненциального распределения и -й степени экспоненциальногоκбИкспе-бИксп+1/(п+1)κΞзнак равноИкс1/(п+1)
κбξпп+1е-бξ/(п+1)1п+1ξ-пп+1знак равноκбп+1е-бξ/(п+1)
ΞE(b/(p+1))κg(x)Е(a)1/(п+1)Е(б/(п+1))распределение, по модулю отсутствующей мультипликативной константы для учета весов: И легко смоделировать как смесь.2
е(Икс)грамм(Икс)знак равно2(12aе-aИкс+12бИкспе-бИксп+1/(п+1))
грамм

Таким образом, R-рендеринг алгоритма принятия-отклонения

simuF <- function(a,b,p){
  reepeat=TRUE
  while (reepeat){
   if (runif(1)<.5) x=rexp(1,a) else
      x=rexp(1,b/(p+1))^(1/(p+1))
   reepeat=(runif(1)>(a+b*x^p)*exp(-a*x-b*x^(p+1)/(p+1))/
      (a*exp(-a*x)+b*x^p*exp(-b*x^(p+1)/(p+1))))}
  return(x)}

и для n-образца:

simuF <- function(n,a,b,p){
  sampl=NULL
  while (length(sampl)<n){
   x=u=sample(0:1,n,rep=TRUE)
   x[u==0]=rexp(sum(u==0),b/(p+1))^(1/(p+1))
   x[u==1]=rexp(sum(u==1),a)
   sampl=c(sampl,x[runif(n)<(a+b*x^p)*exp(-a*x-b*x^(p+1)/(p+1))/
      (a*exp(-a*x)+b*x^p*exp(-b*x^(p+1)/(p+1)))])
   }
  return(sampl[1:n])}

Вот иллюстрация для a = 1, b = 2, p = 3:

введите описание изображения здесь

Сиань
источник