Znalezienie sposobu na symulację liczb losowych dla tego rozkładu

20

Usiłuję napisać program w języku R, który symuluje pseudolosowe liczby z rozkładu za pomocą funkcji rozkładu skumulowanego:

F(x)=1exp(axbp+1xp+1),x0

gdzieza,b>0,p(0,1)

Próbowałem próbkowania z transformacją odwrotną, ale odwrotność nie wydaje się być analitycznie możliwa do rozwiązania. Byłbym zadowolony, gdybyś mógł zasugerować rozwiązanie tego problemu

Sebastian
źródło
1
Za mało czasu na pełną odpowiedź, ale alternatywnie możesz sprawdzić algorytmy próbkowania ważności.
Chuse
1
nie jest to ćwiczenie z podręcznika, ograniczyłem się tylko do tego, ponieważ jest to rozsądne założenie dla moich danych
Sebastian,
6
Jestem zatem zaskoczony „cudowną” normalizacją przez która zmienia rozkład w doskonałą potęgę wykładniczą, ale zdarzają się cuda (z małym prawdopodobieństwem). (p+1)-1
Xi'an

Odpowiedzi:

49

Istnieje proste (i jeśli mogę dodać, eleganckie) rozwiązanie tego ćwiczenia: ponieważ 1-fa(x) wygląda jak iloczyn dwóch rozkładów przeżycia:

(1-fa(x))=exp{-zax-bp+1xp+1}=exp{-zax}1-fa1(x)exp{-bp+1xp+1}1-fa2)(x)
rozkładfajest rozkładem
X=min{X1,X2)}X1fa1,X2)fa2)
W tym przypadkufa1 jestrozkłademwykładniczymmi(za) afa2) jest1/(p+1) -tą potęgąrozkładuwykładniczegomi(b/(p+1)) .

Powiązany kod R jest tak prosty, jak to tylko możliwe

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

i jest zdecydowanie szybszy niż odwrotne pdf i rozdzielczości akceptowania-odrzucania:

> 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 

z zaskakująco idealnym dopasowaniem:

wprowadź opis zdjęcia tutaj

Xi'an
źródło
5
naprawdę fajne rozwiązanie!
Sebastian,
14

Zawsze możesz rozwiązać liczbowo odwrotną transformację.

Poniżej przeprowadzam bardzo proste wyszukiwanie bisekcji. Dla danego prawdopodobieństwa wejściowego q (używam q ponieważ masz już p we wzorze), zaczynam od xL.=0 i xR=1 . Następnie podwajam xRfa(xR)>q . Wreszcie iteracyjnie dzielę przedział [xL.,xR] dopóki jego długość nie będzie mniejsza niż ϵ a jego punkt środkowy xM. spełni fa(xM.)q .

ECDF pasuje do twojego fa wystarczająco dobrze dla moich wyborów za i b , i jest dość szybki. Prawdopodobnie można to przyspieszyć, stosując optymalizację typu Newtona zamiast prostego wyszukiwania bisekcji.

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

S. Kolassa - Przywróć Monikę
źródło
10

Jest nieco skomplikowane, jeśli bezpośrednie rozwiązanie przez akceptację-odrzucenie. Po pierwsze, proste różnicowanie pokazuje, że pdf rozkładu jest

f(x)=(a+bxp)exp{-zax-bp+1xp+1}
f(x)=aeaxebxp+1/(p+1)1+bxpebxp+1/(p+1)eax1
f(x)g(x)=aeax+bxpebxp+1/(p+1)
gξ=xp+1x=ξ1/(p+1)
dxdξ=1p+1ξ1p+11=1p+1ξpp+1
Xκbxpebxp+1/(p+1)κ is the normalising constant, then Ξ=X1/(p+1) has the density
κbξpp+1ebξ/(p+1)1p+1ξpp+1=κbp+1ebξ/(p+1)
which means that (i) Ξ is distributed as an Exponential E(b/(p+1)) variate and (ii) the constant κ is equal to one. Therefore, g(x) ends up being equal to the equally weighted mixture of an Exponential E(a) distribution and the 1/(p+1)-th power of an Exponential E(b/(p+1)) distribution, modulo a missing multiplicative constant of 2 to account for the weights:
fa(x)sol(x)=2)(12)zami-zax+12)bxpmi-bxp+1/(p+1))
I sol można łatwo symulować jako mieszaninę.

Renderowanie R algorytmu akceptacji-odrzucenia jest zatem

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)}

i dla próbki 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])}

Oto ilustracja dla a = 1, b = 2, p = 3:

wprowadź opis zdjęcia tutaj

Xi'an
źródło