Menghasilkan sampel acak dari distribusi khusus

Saya mencoba untuk menghasilkan sampel acak dari pdf kustom menggunakan R. Pdf saya adalah:

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$

Saya menghasilkan sampel yang seragam dan kemudian mencoba mengubahnya menjadi distribusi khusus saya. Saya melakukan ini dengan menemukan cdf dari distribusi saya ( $F_{X}(x)$ ) dan mengaturnya ke sampel seragam ( ) dan menyelesaikannya untuk . $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})$

Untuk menghasilkan sampel acak dengan distribusi di atas, dapatkan sampel seragam dan selesaikan untuk in $u \in[0,1]$ $x$

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

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

Saya menerapkannya Rdan saya tidak mendapatkan distribusi yang diharapkan. Adakah yang bisa menunjukkan kekurangan dalam pemahaman saya?

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
sumber

Pasti ada kesalahan pengkodean. Saya tidak menggunakan R, jadi saya tidak bisa mengatakan apa kesalahan sebenarnya - tapi saya baru saja membuat kode solusi Anda (berhati-hati untuk mengambil akar tengah dari polinomial kubik, yang selalu terletak antara 0 dan 1), dan Saya mendapatkan persetujuan yang baik antara sampel dan distribusi yang diharapkan. Mungkinkah ada masalah dengan pencari root Anda? Apa yang salah dengan sampel yang Anda dapatkan?

jpillow

Saya mencoba kode Anda (yang tidak sangat efisien, omong-omong) dan mendapatkan distribusi yang diharapkan.

Aniko

@jpillow dan @Aniko Kesalahan saya. Ketika saya menggunakan nsamples <- 1e6itu cocok.

Anand

@ Anand Salah satu cara adalah dengan mengamati bahwa , memungkinkan perhitungan langsung dalam hal .

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

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

x

$x$

u

$u$

whuber

@Anand en.wikipedia.org/wiki/…

whuber

Jawaban:

Sepertinya Anda tahu bahwa kode Anda berfungsi, tetapi @Aniko menunjukkan bahwa Anda dapat meningkatkan efisiensinya. Kenaikan kecepatan terbesar Anda mungkin berasal dari pra-alokasi memori zagar Anda tidak menanamnya di dalam satu lingkaran. Sesuatu seperti z <- rep(NA, nsamples)harus melakukan triknya. Anda dapat memperoleh kecepatan kecil dengan menggunakan vapply()(yang menentukan tipe variabel yang dikembalikan) alih-alih loop eksplisit (ada pertanyaan SO yang bagus tentang keluarga yang berlaku).

> 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

Dan Anda tidak perlu ;di akhir setiap baris (apakah Anda seorang konversi MATLAB?).

Richard Herron
sumber

Terima kasih atas jawaban terperinci Anda dan untuk menunjukkannya vapply. Saya telah mengkode dalam C/C++waktu yang sangat lama dan itulah alasan ;kesengsaraan!

Anand

uniroot

10^{7}

$10^7$