Pendahuluan

Metode Polar-Marsaglia merupakan salah satu metode untuk menghasilkan sepasang peubah acak normal baku yang independen. Metode ini merupakan variasi dari metode Box-Muller. Perbedaannya terletak pada penggunaan komputasi:
- Box-Muller menggunakan fungsi trigonometri seperti sinus dan cosinus,
- Sedangkan Polar-Marsaglia menghindari komputasi trigonometri dengan menggunakan transformasi langsung melalui perhitungan kuadrat dari variabel acak.

Persamaan transformasi pada metode Polar-Marsaglia adalah sebagai berikut:

Persamaan Box-Muller

\(N_1 = \sqrt{-2 \ln U_1} \cos(2\pi U_2)\)

\(N_2 = \sqrt{-2 \ln U_1} \sin(2\pi U_2)\)

Bentuk alternatif (dengan R^2)

\(N_1 = \sqrt{-2 \ln R^2}\; V_1 \; (V_1^2+V_2^2)^{-1/2}\)

\(N_2 = \sqrt{-2 \ln R^2}\; V_2 \; (V_1^2+V_2^2)^{-1/2}\)

Bentuk yang disederhanakan

\(N_1 = \sqrt{-2 \ln (V_1^2+V_2^2)}\; V_1 \; (V_1^2+V_2^2)^{-1/2}\)

\(N_2 = \sqrt{-2 \ln (V_1^2+V_2^2)}\; V_2 \; (V_1^2+V_2^2)^{-1/2}\)

Bentuk akhir

\(N_1 = V_1 \sqrt{\frac{-2 \ln W}{W}}, \qquad N_2 = V_2 \sqrt{\frac{-2 \ln W}{W}},\) dengan \(W = V_1^2 + V_2^2.\)

dengan:

- \(( V_1 )\) dan \(( V_2 )\) merupakan transformasi dari variabel acak uniform,

- \(( W = V_1^2 + V_2^2 )\),

- Nilai \(( W )\) harus kurang dari 1 agar transformasi valid.

Implementasi di R

Inisialisasi Variabel

Pertama, kita inisialisasi beberapa variabel dasar dan menghasilkan bilangan acak uniform.

# Jumlah bilangan acak yang dihasilkan
i <- 1000

# Mengatur seed untuk reproduktifitas
set.seed(123)

# Menghasilkan dua vektor bilangan acak uniform
U1 <- runif(i)
U2 <- runif(i)

# Transformasi ke rentang [-1, 1]
V1 <- 2 * U1 - 1
V2 <- 2 * U2 - 1

# Menghitung W = V1^2 + V2^2
W <- V1^2 + V2^2

Seleksi Nilai yang Valid

Hanya nilai-nilai dengan \(W < 1\) yang valid untuk proses selanjutnya. Jika jumlah nilai yang valid kurang dari jumlah yang diinginkan, maka dilakukan pengulangan untuk mengisi kekurangan.

# Menyaring nilai dengan W < 1
W1 <- W[W < 1]
V1 <- V1[W < 1]
V2 <- V2[W < 1]

# Menghitung jumlah nilai yang kurang
i_1 <- i - length(W1)

# Lakukan perulangan sampai jumlah nilai yang valid mencapai jumlah i
while (i_1 > 0) {
  U_1 <- runif(i_1)
  U_2 <- runif(i_1)
  V_1 <- 2 * U_1 - 1
  V_2 <- 2 * U_2 - 1
  W_1 <- V_1^2 + V_2^2
  
  # Menyaring nilai baru yang valid
  valid <- W_1 < 1
  W1 <- c(W1, W_1[valid])
  V1 <- c(V1, V_1[valid])
  V2 <- c(V2, V_2[valid])
  
  # Update jumlah nilai yang kurang
  i_1 <- i - length(W1)
}

Transformasi ke Distribusi Normal

Setelah mendapatkan nilai yang valid, transformasi dilakukan untuk mendapatkan sepasang bilangan acak yang berdistribusi normal standar.

# Transformasi ke bilangan acak normal
N1 <- V1 * sqrt(-2 * log(W1) / W1)
N2 <- V2 * sqrt(-2 * log(W1) / W1)

# Menampilkan statistik dasar dari N1 dan N2
hasil <- data.frame(
  parameter = c("miu N1", "sigma2 N1", "miu N2", "sigma2 N2"),
  Nilai = c(mean(N1), var(N1), mean(N2), var(N2))
)
hasil
##   parameter        Nilai
## 1    miu N1  0.008316551
## 2 sigma2 N1  1.062962164
## 3    miu N2 -0.012600728
## 4 sigma2 N2  0.979850242

Visualisasi

Berikut adalah histogram dari bilangan acak normal yang dihasilkan, yang diharapkan mendekati distribusi Normal \(\mathcal{N}(0,1)\).

par(mfrow = c(1, 2))
hist(N1, col = "Thistle", main = "Histogram of N1 ~ Normal (0,1)")
hist(N2, col = "Thistle", main = "Histogram of N2 ~ Normal (0,1)")