Generowanie losowych próbek z niestandardowej dystrybucji

16

Próbuję wygenerować losowe próbki z niestandardowego pliku PDF przy użyciu R. Mój pdf to:

fX(x)=32(1x2),0x1

Wygenerowałem jednolite próbki, a następnie próbowałem przekształcić je w moją niestandardową dystrybucję. Zrobiłem to, znajdując plik cdf mojej dystrybucji ( FX(x) ) i ustawiając go na jednolitą próbkę ( ) i rozwiązując dla .ux

FX(x)=Pr[Xx]=0x32(1y2)dy=32(xx33)

Aby wygenerować losową próbkę z powyższym rozkładem, pobierz jednolitą próbkę u[0,1] i rozwiąż dla x in

32(xx33)=u

Zaimplementowałem to w Ri nie otrzymuję oczekiwanej dystrybucji. Czy ktoś może wskazać wadę w moim rozumieniu?

nsamples <- 1000;
x <- runif(nsamples);

f <- function(x, u) { 
  return(3/2*(x-x^3/3) - u);
}

z <- c();
for (i in 1:nsamples) {
  # find the root within (0,1) 
  r <- uniroot(f, c(0,1), tol = 0.0001, u = x[i])$root;
  z <- c(z, r);
}
Anand
źródło
1
To musi być błąd kodowy. Nie używam R, więc nie mogę powiedzieć, jaki dokładnie jest błąd - ale właśnie zakodowałem twoje rozwiązanie (uważając, aby wziąć środkowy pierwiastek wielomianu sześciennego, który zawsze leży między 0 a 1), i Otrzymuję dobrą zgodność między próbkami a oczekiwanym rozkładem. Czy to może być problem z twoją root rootką? Co jest nie tak z próbkami, które otrzymujesz?
jpillow
Wypróbowałem twój kod (który, nawiasem mówiąc, nie jest zbyt wydajny) i dostaję oczekiwaną dystrybucję.
Aniko,
@jpillow i @Aniko Mój błąd. Kiedy użyłem nsamples <- 1e6, było to dobre dopasowanie.
Anand
2
@Anand Jednym ze sposobów jest zaobserwowanie, że , pozwalając na bezpośrednie obliczenie pod względem . x=2sin(arcsin(u)/3)xu
whuber
1
@Anand en.wikipedia.org/wiki/…
whuber

Odpowiedzi:

11

Wygląda na to, że zorientowałeś się, że Twój kod działa, ale @Aniko wskazało, że możesz poprawić jego wydajność. Twój największy wzrost prędkości prawdopodobnie pochodziłby z alokacji pamięci dla z, abyś nie rozwijał jej w pętli. Coś takiego z <- rep(NA, nsamples)powinno załatwić sprawę. Możesz uzyskać niewielki wzrost prędkości przy użyciu vapply()(który określa zwracany typ zmiennej) zamiast wyraźnej pętli (istnieje świetne pytanie SO dotyczące rodziny zastosuj).

> nsamples <- 1E5
> x <- runif(nsamples)
> f <- function(x, u) 1.5 * (x - (x^3) / 3) - u
> z <- c()
> 
> # original version
> system.time({
+ for (i in 1:nsamples) {
+   # find the root within (0,1) 
+   r <- uniroot(f, c(0,1), tol = 0.0001, u = x[i])$root
+   z <- c(z, r)
+ }
+ })
   user  system elapsed 
  49.88    0.00   50.54 
> 
> # original version with pre-allocation
> z.pre <- rep(NA, nsamples)
> system.time({
+ for (i in 1:nsamples) {
+   # find the root within (0,1) 
+   z.pre[i] <- uniroot(f, c(0,1), tol = 0.0001, u = x[i])$root
+   }
+ })
   user  system elapsed 
   7.55    0.01    7.78 
> 
> 
> 
> # my version with sapply
> my.uniroot <- function(x) uniroot(f, c(0, 1), tol = 0.0001, u = x)$root
> system.time({
+   r <- vapply(x, my.uniroot, numeric(1))
+ })
   user  system elapsed 
   6.61    0.02    6.74 
> 
> # same results
> head(z)
[1] 0.7803198 0.2860108 0.5153724 0.2479611 0.3451658 0.4682738
> head(z.pre)
[1] 0.7803198 0.2860108 0.5153724 0.2479611 0.3451658 0.4682738
> head(r)
[1] 0.7803198 0.2860108 0.5153724 0.2479611 0.3451658 0.4682738

I nie potrzebujesz ;znaku na końcu każdej linii (czy jesteś konwerterem MATLAB?).

Richard Herron
źródło
Dziękujemy za szczegółową odpowiedź i wskazanie vapply. Wchodzę w C/C++to od bardzo dawna i to jest przyczyną ;dolegliwości!
Anand
1
uniroot107