Generowanie losowych próbek z niestandardowej dystrybucji

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

f_{X} (x) = \frac{3}{2} (1 - x^{2}), 0 \leq x \leq 1

$f_{X}(x) = \frac{3}{2} (1-x^2), 0 \le x \le 1$

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 ( $F_{X}(x)$ ) i ustawiając go na jednolitą próbkę ( ) i rozwiązując dla . $u$ $x$

F_{X} (x) = Pr [X \leq x] = \int_{0}^{x} \frac{3}{2} (1 - y^{2}) d y = \frac{3}{2} (x - \frac{x^{3}}{3})

$F_{X}(x) = \Pr[X \le x] = \int_{0}^{x} \frac{3}{2} (1-y^2) dy = \frac{3}{2} (x - \frac{x^3}{3})$

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

\frac{3}{2} (x - \frac{x^{3}}{3}) = u

$\frac{3}{2} (x - \frac{x^3}{3}) = 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);
}

r sampling uniform

— Anand
źródło

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

@Anand Jednym ze sposobów jest zaobserwowanie, że , pozwalając na bezpośrednie obliczenie pod względem .

x = 2 \sin (\arcsin (u) / 3)

$x = 2\sin(\arcsin(u)/3)$

x

$x$

u

$u$

— whuber

@Anand en.wikipedia.org/wiki/…

— whuber

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

uniroot

10^{7}

$10^7$