Generación de muestras aleatorias de una distribución personalizada

Question

Estoy tratando de generar muestras aleatorias a partir de un pdf personalizados mediante R. Mi pdf es: $$f_{X}(x) = \frac{3}{2} (1-x^2), 0 \le x \le 1$$

I generated uniform samples and then tried to transform it to my custom distribution. I did this by finding the cdf of my distribution ($F_{X}(x)$) and setting it to the uniform sample ($u$) and solving for $x$.

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

To generate a random sample with the above distribution, get a uniform sample $u \en[0,1]$ and solve for $x$ in $$\frac{3}{2} (x - \frac{x^3}{3}) = u $$

He implementado en R y no me da la distribución esperada. Puede alguien señalar el error en mi entendimiento?

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

Answer 1

Parece que descubriste que el código funciona, pero @Aniko señaló que se podría mejorar su eficiencia. Su mayor ganancia de velocidad probablemente provienen de pre-asignación de memoria para z , de modo que usted no está creciendo dentro de un bucle. Algo como z <- rep(NA, nsamples) debe hacer el truco. Usted puede obtener una pequeña ganancia de velocidad de uso vapply() (en el que se especifica que devuelve el tipo de variable) en lugar de un explícito de bucle (hay una gran ASÍ, pregunta en el aplique de la familia).

> 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

Y usted no necesita la ; al final de cada línea (son una MATLAB convertir?).