Suche nach einer Möglichkeit, Zufallszahlen für diese Verteilung zu simulieren

20

Ich versuche, ein Programm in R zu schreiben, das Pseudozufallszahlen aus einer Verteilung mit der kumulativen Verteilungsfunktion simuliert:

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

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

Ich habe versucht, die inverse Transformation abzutasten, aber die inverse scheint analytisch nicht lösbar zu sein. Ich würde mich freuen, wenn Sie eine Lösung für dieses Problem vorschlagen könnten

Sebastian
quelle
1
Nicht genügend Zeit für eine vollständige Antwort, aber Sie können alternativ auch die Algorithmen für das wichtige Sampling überprüfen.
Chuse
1
es ist keine Lehrbuchübung, ich habe die Einschränkung nur festgelegt, weil dies eine vernünftige Annahme für meine Daten ist
Sebastian
6
Ich bin dann überrascht über die "wundersame" Normalisierung durch , die die Verteilung in eine perfekte Potenz eines Exponentials verwandelt, aber Wunder geschehen (mit geringer Wahrscheinlichkeit). (p+1)1
Xi'an

Antworten:

49

Es gibt eine einfache (und wenn ich hinzufügen darf, elegante) Lösung für diese Aufgabe: da wie ein Produkt zweier Überlebensverteilungen erscheint: die Verteilung ist die Verteilung von In diesem Fall ist F_1 die Exponential \ mathcal {E} (a) -Verteilung und F_2 ist die 1 / (p + 1) -te Potenz einer Exponential \ mathcal {E} (b / (p + 1)) - Verteilung.1F(x)

(1F(x))=exp{axbp+1xp+1}=exp{ax}1F1(x)exp{bp+1xp+1}1F2(x)
F
X=min{X1,X2}X1F1,X2F2
F1E(a)F21/(p+1)E(b/(p+1))

Der zugehörige R-Code ist so einfach wie es nur geht

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

und es ist definitiv viel schneller als das Inverse PDF und das Akzeptieren-Ablehnen von Auflösungen:

> 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 

mit einer nicht überraschend perfekten Passform:

Bildbeschreibung hier eingeben

Xi'an
quelle
5
wirklich coole lösung!
Sebastian
14

Sie können die inverse Transformation immer numerisch lösen.

Im Folgenden führe ich eine sehr einfache Halbierungssuche durch. Für eine gegebene Eingabewahrscheinlichkeit (ich verwende da Sie bereits ein in Ihrer Formel haben) beginne ich mit und . Dann ich bis . Schließlich halbiere ich iterativ das Intervall bis seine Länge kürzer als und sein Mittelpunkt erfüllt .qqpxL=0xR=1xRF(xR)>q[xL,xR]ϵxMF(xM)q

Das ECDF passt gut genug zu Ihrem für meine Auswahl von und , und es ist ziemlich schnell. Sie könnten dies wahrscheinlich beschleunigen, indem Sie anstelle der einfachen Halbierungssuche eine Optimierung nach Newton verwenden.Fab

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 - Setzen Sie Monica wieder ein
quelle
10

Es gibt eine etwas verworrene, wenn direkte Auflösung durch Akzeptieren-Ablehnen. Zunächst zeigt eine einfache Unterscheidung, dass das PDF der Distribution Zweitens, da Wir haben die obere Schranke drittes nehme man unter Berücksichtigung des zweiten Terms in die Änderung der Variablen , dh . Dann ist ist der Jacobi der Variablenänderung. Wenn

f(x)=(a+bxp)exp{axbp+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
Xhat eine Dichte der Form wobei die Normalisierungskonstante ist, dann ist hat die Dichte was bedeutet, dass (i) ist verteilt als Exponential variate und (ii) die Konstante ist gleich eins. Daher ist gleich der gleichgewichteten Mischung einer Exponential -Verteilung und der -ten Potenz einer Exponentialκbxpebxp+1/(p+1)κΞ=X1/(p+1)
κbξpp+1ebξ/(p+1)1p+1ξpp+1=κbp+1ebξ/(p+1)
ΞE(b/(p+1))κg(x)E(a)1/(p+1)E(b/(p+1))Verteilung, modulo eine fehlende multiplikative Konstante von , um die Gewichte zu berücksichtigen: Und ist einfach als Mischung zu simulieren.2
f(x)g(x)=2(12aeax+12bxpebxp+1/(p+1))
g

Ein R-Rendering des Accept-Reject-Algorithmus liegt somit vor

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

und für eine n-Probe:

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

Hier ist eine Illustration für a = 1, b = 2, p = 3:

Bildbeschreibung hier eingeben

Xi'an
quelle