Génération d'échantillons aléatoires à partir d'une distribution personnalisée

J'essaie de générer des échantillons aléatoires à partir d'un pdf personnalisé à l'aide de R. Mon pdf est:

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$

J'ai généré des échantillons uniformes, puis j'ai essayé de le transformer en ma distribution personnalisée. J'ai fait cela en trouvant le cdf de ma distribution ( ) et en le réglant sur l'échantillon uniforme ( ) et en résolvant pour . $F_{X}(x)$ $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})$

Pour générer un échantillon aléatoire avec la distribution ci-dessus, obtenez un échantillon uniforme et résolvez pour en $u \in[0,1]$ $x$

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

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

Je l'ai implémenté dans Ret je n'obtiens pas la distribution attendue. Quelqu'un peut-il signaler la faille dans ma compréhension?

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
la source

Doit être une erreur de codage. Je n'utilise pas R, donc je ne peux pas dire exactement quelle est l'erreur - mais je viens de coder votre solution (en prenant soin de prendre la racine centrale du polynôme cubique, qui se situe toujours entre 0 et 1), et J'obtiens un bon accord entre les échantillons et la distribution attendue. Serait-ce un problème avec votre root finder? Quel est le problème avec les échantillons que vous obtenez?

jpillow

J'ai essayé votre code (qui n'est pas très efficace d'ailleurs) et j'obtiens la distribution attendue.

Aniko

@jpillow et @Aniko Mon erreur. Quand je l'ai utilisé, nsamples <- 1e6c'était un bon match.

Anand

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

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

x

$x$

u

$u$

@Anand en.wikipedia.org/wiki/…

whuber

Réponses:

Il semble que vous ayez compris que votre code fonctionne, mais @Aniko a souligné que vous pouviez améliorer son efficacité. Votre gain de vitesse le plus important proviendrait probablement de la pré-allocation de mémoire pour zque vous ne la développiez pas dans une boucle. Quelque chose comme z <- rep(NA, nsamples)ça devrait faire l'affaire. Vous pouvez obtenir un petit gain de vitesse en utilisant vapply()(qui spécifie le type de variable renvoyé) au lieu d'une boucle explicite (il y a une grande question SO sur la famille apply).

> 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

Et vous n'avez pas besoin du ;à la fin de chaque ligne (êtes-vous un converti MATLAB?).

Richard Herron
la source

Merci pour votre réponse détaillée et pour avoir souligné vapply. J'encode C/C++depuis très longtemps et c'est la raison de l' ;affliction!

Anand

+1 Il semble que le remplacement unirootpar une solution directe accélère les choses de trois ordres de grandeur supplémentaires: vous ne devriez avoir aucun problème à vous déplacer

10^{7}

$10^7$ varie par seconde.

whuber