Usiłuję napisać program w języku R, który symuluje pseudolosowe liczby z rozkładu za pomocą funkcji rozkładu skumulowanego:
gdzie
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
r
random-generation
Sebastian
źródło
źródło
Odpowiedzi:
Istnieje proste (i jeśli mogę dodać, eleganckie) rozwiązanie tego ćwiczenia: ponieważ1−F(x) wygląda jak iloczyn dwóch rozkładów przeżycia:
( 1 - F( x ) ) = exp{ - a x - bp + 1xp + 1} = exp{ - a x }1 - F1( x )exp{ - bp + 1xp + 1}1 - F2)( x )
rozkładfa jest rozkładem
X= min { X1, X2)}X1∼ F.1, X2)∼ F.2)
W tym przypadkufa1 jestrozkłademwykładniczymmi( ) 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
i jest zdecydowanie szybszy niż odwrotne pdf i rozdzielczości akceptowania-odrzucania:
z zaskakująco idealnym dopasowaniem:
źródło
Zawsze możesz rozwiązać liczbowo odwrotną transformację.
Poniżej przeprowadzam bardzo proste wyszukiwanie bisekcji. Dla danego prawdopodobieństwa wejściowegoq (używam q ponieważ masz już p we wzorze), zaczynam od xL.= 0 i xR= 1 . Następnie podwajam xR aż fa( 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 twojegofa 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.
źródło
Jest nieco skomplikowane, jeśli bezpośrednie rozwiązanie przez akceptację-odrzucenie. Po pierwsze, proste różnicowanie pokazuje, że pdf rozkładu jestf(x)=(a+bxp)exp{−ax−bp+1xp+1} f(x)=ae−axe−bxp+1/(p+1)≤1+bxpe−bxp+1/(p+1)e−ax≤1 f(x)≤g(x)=ae−ax+bxpe−bxp+1/(p+1) g ξ=xp+1 x=ξ1/(p+1) dxdξ=1p+1ξ1p+1−1=1p+1ξ−pp+1 X κbxpe−bxp+1/(p+1) κ is the normalising constant, then Ξ=X1/(p+1) has the density
κbξpp+1e−bξ/(p+1)1p+1ξ−pp+1=κbp+1e−bξ/(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:
f(x)≤g( x ) = 2 ( 12)e- a x+ 12)b xpmi- b xp + 1/ (p+1))
I sol można łatwo symulować jako mieszaninę.
Renderowanie R algorytmu akceptacji-odrzucenia jest zatem
i dla próbki n:
Oto ilustracja dla a = 1, b = 2, p = 3:
źródło