Zufallsstichprobe aus inversem Cdf ohne geschlossene Form generieren

Ich arbeite an einer bestimmten Distribution, deren inverses cdf nicht in geschlossener Form existiert. Das cdf der Distribution ist gegeben durch

F (x; d, m, p, α, β) = \frac{1 - (1 + x^{m})^{- d} \exp (- β x^{α})}{1 - p (1 + x^{m})^{- d} \exp (- β x^{α})}

$F(x; d, m, p, \alpha, \beta) = \frac{1-(1+x^m)^{-d} \exp(-\beta x^\alpha)}{1-p(1+x^m)^{-d} \exp(-\beta x^\alpha)}$

für streng positive und . $m, d, \alpha, \beta$ $0\lt p \lt 1$

Mein Problem ist, dass ich neu im RPaket bin und eine Zufallsstichprobe aus der Verteilung mit generieren muss R.

r distributions random-generation Soma
quelle

NB: Diese Formel für gibt nicht immer ein cdf: Sie müssen es auf die Domäne damit es definiert und gültig ist.

F

$F$

x \geq 0

$x\ge 0$

whuber

Entschuldigung, ich habe vergessen, es auszudrücken. Es ist cdf mit x≥0. Danke @Whuber.

Soma

Um zu invertieren, müssen Sie numerisch lösen, da dies eine transzendentale Gleichung ist. Vielleicht ist hier ein cleveres Akzeptanz- / Ablehnungsschema angebracht?

F

$F$

y = \exp (- β x^{α}) / (1 + x^{m})

$y=\exp(-\beta x^\alpha)/(1+x^m)$

StijnDeVuyst

Hast du diesen Thread

user603

Das direkte Lösen der Gleichung kann für die höchsten Quantile problematisch und unsicher sein . Löse stattdessen für und nimm . Mit einer kleinen Anzahl von Newton-Raphson-Iterationen erhalten Sie das Ergebnis mit hoher Präzision. (Dies ist so einfach, dass Sie es sogar in einer Tabelle implementieren können.)

q

$q$

\log (1 - q) = \log (1 - F (x; d, m, p, 1, 1)) = \log (\frac{1 - p}{- p + e^{z} {(z^{m} + 1)}^{d}})

$\log(1-q) =\log(1-F(x;d,m,p,1,1))= \log\left(\frac{1-p}{-p+e^z \left(z^m+1\right)^d}\right)$

z \geq 0

$z \ge 0$

x = (z / β)^{1 / α}

$x = (z/\beta)^{1/\alpha}$

whuber

Antworten:

Hier sind zwei Möglichkeiten, um numerische Näherungen an die Umkehrung des cdf zu berechnen, vorausgesetzt, Sie haben Entscheidungen für m, d, α, β und p getroffen. Beide Methoden erfordern, dass Sie F (x) für ein gegebenes x berechnen können, also ...

m = 1
d = 2
a = 1
b = 2
p = 0.5
F = function(x) (1 - ((1+x^m)^-d) * exp(-b*x^a))  / 
    (1 - (p*(1+x^m)^-d) * exp(-b*x^a))

Methode 1

Um InvF (a) zu berechnen, lösen Sie die Gleichung F (x) = a

InvF1 = function(a) uniroot(function(x) F(x) - a, c(0,10))$root
InvF1(0.5)
[1] 0.1038906
F(InvF1(0.5))
[1] 0.4999983

Methode 2

Bewerten Sie y = F (x) für einen Bereich von x und passen Sie dann eine Kurve als Funktion von y an x an.

x = c(seq(0,3, 0.001), seq(3.1,10,0.1))
y = F(x)
InvF2 = approxfun(y, x)

InvF2(0.5)
[1] 0.1038916
F(InvF2(0.5))
[1] 0.5000011

Sie können die Genauigkeit erhöhen, InvF2indem Sie eine dichtere Abtastung von x verwenden, insbesondere für kleine Werte von x.

G5W
quelle