Pembangkitan Bilangan Acak
Hazelita Dwi Rahmasari
2024-09-19
\[f(x) = \frac{1}{\beta} \exp\left(-\frac{(x-\mu)}{\beta} - e^{-\frac{(x-\mu)}{\beta}}\right), -\infty < x < \infty\]
#mendefinisikan nilai (mu, beta)
set.seed(123)
mu = 0 # Nilai rata-rata (mu)
beta = 1 # Skala distribusi (beta)# Fungsi untuk inverse CDF dari distribusi Gumbel
invers.fx = function(u) {
mu - beta * log(-log(u))
}
# Fungsi pembangkit bilangan acak menggunakan inverse transform
rand.fx = function(n, fxu) {
u = runif(n) # Bangkitkan n bilangan acak dari uniform(0,1)
x = fxu(u) # Terapkan fungsi fyu (di sini fungsi inverse CDF)
return(x) # Kembalikan nilai y
}# Bangkitkan 100 bilangan acak dari distribusi Gumbel
set.seed(123)
n = 100
fx = rand.fx(n, invers.fx)
fx## [1] -0.22014933 1.43603082 0.11194151 2.08416975 2.79069703 -1.12778331
## [7] 0.44869746 2.17314163 0.51880524 0.24345433 3.12070507 0.23429832
## [13] 0.94355538 0.58427564 -0.82143388 2.24852302 -0.33794769 -1.15331169
## [19] -0.10883953 3.06693159 2.14514065 1.00236889 0.80856809 5.15912879
## [25] 0.86264791 1.06548060 0.49645493 0.65270242 -0.21573701 -0.65052663
## [31] 3.27871331 2.27487967 0.99413860 1.47479605 -1.30953553 0.30303740
## [37] 1.28566806 -0.42565332 -0.13552240 -0.38023769 -0.66593531 0.12718526
## [43] 0.12493371 0.00262589 -0.63177859 -0.68040512 -0.37608469 0.26964541
## [49] -0.28093071 1.87501936 -1.12583533 0.20334970 1.49393125 -0.74410657
## [55] 0.54796154 -0.45571634 -0.72241017 1.26131494 2.19929879 0.01789631
## [61] 0.89699061 -0.85677724 0.04375149 -0.25714141 1.58470172 0.22088340
## [67] 1.55759723 1.57129777 1.46862967 0.19678994 1.26679573 0.76943979
## [73] 1.07226221 -1.99851909 0.29601775 -0.41448266 0.03245404 0.71383168
## [79] -0.04372822 -0.78709543 -0.34511181 0.90786658 0.13568316 1.43544805
## [85] -0.82169047 0.18313487 4.18927210 2.17938157 2.11603039 -0.55541755
## [91] -0.71043863 0.85326412 -0.06627471 0.86645475 -0.12950832 -0.51459311
## [97] 1.40435950 -0.86237440 0.27194091 0.39988518# Buat histogram dari bilangan acak yang dihasilkan
hist(fx, prob=TRUE, col = "steelblue", xlab = "x",
main=expression(f(y)==(1/beta) * exp(-(x - mu)/beta - exp(-(x - mu)/beta))))
# Tambahkan FKP (PDF) ke histogram
curve((1/beta) * exp(-(x - mu)/beta - exp(-(x - mu)/beta)),
col="red", lwd=2, add=TRUE)
#Acceptance-Rejection Method
# Fungsi untuk Acceptance-Rejection Method
x2 = function(x0, x1, n, f) {
xx = seq(x0, x1, len = 1001)
mf = max(f(xx)) # Nilai maksimum dari PDF Gumbel dalam interval [x0, x1]
k = 0
j = 0
hasil = numeric(n)
# Loop hingga jumlah n tercapai
while (k < n) {
x = runif(1, x0, x1) # Bilangan acak dari distribusi uniform di [x0, x1]
j = j + 1
y1 = runif(1, 0, mf) # Bilangan acak dari distribusi uniform di [0, mf]
y2 = f(x) # Evaluasi PDF Gumbel di titik u
# Jika nilai fungsi PDF lebih besar dari y1, terima sampel
if (y2 >= y1) {
k = k + 1
hasil[k] = x
}
}
# Kembalikan hasil dan jumlah iterasi
list(hasil=hasil, j=j)
}set.seed(123)
f = function(x) {(1 / beta) * exp(-(x - mu) / beta - exp(-(x - mu) / beta))}
x1 = x2(-10, 10, 100, f)
x1## $hasil
## [1] 1.02870029 0.88132049 1.21895968 3.30230389 4.20364803 -0.49366852
## [7] -0.66441917 -0.22773933 -0.34195206 1.70966706 0.42271452 1.43870628
## [13] -1.10463996 0.04599127 2.59946107 -0.35914789 -0.44309240 0.84160734
## [19] 0.93652313 1.86091303 0.62140798 0.32889789 -0.95083212 0.52059320
## [25] -0.27176472 0.30143616 -0.43087232 0.94918936 1.92527089 0.47645203
## [31] -0.71966746 1.69873115 2.54292369 0.72107441 -0.70572452 0.37609847
## [37] 0.71082783 0.59428183 0.19685799 -1.73389184 1.33468166 0.53180159
## [43] 0.16741259 2.59490848 -0.48498989 0.05816608 -0.81150657 0.44671532
## [49] -0.41456880 0.79121828 0.92852517 -0.26354253 1.14875166 -0.11415276
## [55] 0.27562278 1.38705377 -0.58636333 -0.23009035 -0.97662187 2.25446937
## [61] 1.00518093 0.08879622 -1.51095810 0.15365726 0.51948548 1.20561075
## [67] 2.90855031 0.16316169 2.46522767 0.96561148 1.44297359 -0.75802284
## [73] -0.63135763 1.71131153 -1.04859718 -0.86226489 1.10627184 1.01845458
## [79] -0.40911844 0.63366654 0.57530792 0.20822603 -0.81518366 -0.21817003
## [85] -0.12767623 1.41426387 1.35452332 1.25033617 0.96086652 -0.22652930
## [91] 0.25149943 -1.23450925 0.72719862 -0.72852472 2.27681887 1.79051301
## [97] -0.74602317 0.01343398 -0.44702756 -0.87360878
##
## $j
## [1] 752# Buat histogram dari nilai yang dihasilkan
hist(x1[[1]], prob=TRUE,
main=expression((1 / beta) * exp(-(x - mu) / beta - exp(-(x - mu) / beta))),
xlab="x", col="yellow")
# Buat kurva PDF Gumbel di atas histogram
tx = seq(-10, 10, 0.01)
lines(tx, f(tx), col=2, lwd=2) # Tambahkan PDF ke histogram
par(mfrow=(c(1,2)))
hist(fx, prob=TRUE, col = "steelblue", xlab="x", main="Inverse
Transform Method")
curve((1/beta) * exp(-(x - mu)/beta - exp(-(x - mu)/beta)),
col="red", lwd=2, add=TRUE)
#############################################
hist(x1[[1]], prob=TRUE,
main="Acceptance-Rejection Method",
xlab="x", col="yellow", border="blue")
tx = seq(-10, 10, 0.01)
lines(tx, f(tx), col=2, lwd=2) 
\[ f(y) = \begin{cases} y & \text{untuk } 0 \leq y \leq 1 \\ 1 & \text{untuk } 1 < y \leq 1.5 \\ 0 & \text{lainnya} \end{cases} \]
set.seed(123)
# Fungsi Invers Transform Method (ITM)
fy = function(n, fyu) {
u = runif(n) # Bilangan acak dari distribusi uniform(0,1)
y = fyu(u) # Terapkan fungsi invers CDF
return(y) # Kembalikan nilai y
}
# Invers CDF untuk fungsi f(y)
invers_b = function(u) {
ifelse(u <= 0.5, sqrt(2 * u), 1 + (u - 0.5))
}
# Menghasilkan 100 bilangan acak menggunakan Invers Transform Method
bilangan_acak = fy(100, invers_b)
bilangan_acak## [1] 0.75838977 1.28830514 0.90440801 1.38301740 1.44046728 0.30184930
## [7] 1.02810549 1.39241904 1.05143501 0.95563040 1.45683335 0.95219132
## [13] 1.17757064 1.07263340 0.45370625 1.39982497 0.70155219 0.29003287
## [19] 0.80984038 1.45450365 1.38953932 1.19280341 1.14050681 1.49426978
## [25] 1.15570580 1.20853047 1.04406602 1.09414202 0.76047319 0.54242723
## [31] 1.46302423 1.40229905 1.19070528 1.29546742 0.22187242 0.97754383
## [37] 1.25845954 0.65788743 0.79772302 0.68062587 0.53441561 0.91054526
## [43] 0.90964205 0.85888934 0.55216800 0.52688910 0.68269188 0.96536257
## [49] 0.72934579 1.35782772 0.30275788 0.94042551 1.29892485 0.49375958
## [55] 1.06094798 0.64269960 0.50503792 1.25330786 1.39504536 0.86540485
## [61] 1.16511519 0.43552419 0.87632145 0.74078829 1.31464004 0.94711809
## [67] 1.31006435 1.31238951 1.29434232 0.93790371 1.25447516 1.12922113
## [73] 1.21018240 0.03534893 0.97500418 0.66350416 0.87156932 1.11277100
## [79] 0.83880619 0.47145609 0.69802503 1.16805559 0.91394396 1.28819583
## [85] 0.45357391 0.93262291 1.48495698 1.39305111 1.38646906 0.59169697
## [91] 0.51126449 1.15310193 0.82887450 1.15675813 0.80046642 0.61268445
## [97] 1.28229430 0.43265457 0.96620810 1.01150546# Acceptance-Rejection Method (ARM) yang langsung mengembalikan bilangan acak
ARM_simplified = function(x0, x1, n, f) {
xx = seq(x0, x1, len = 1001) # Rentang kandidat
mf = max(f(xx)) # Maksimum f(y)
k = 0 # Counter bilangan acak yang diterima
hasil = numeric(n) # Tempat penyimpanan bilangan acak
while (k < n) {
u = runif(1, x0, x1) # Kandidat dari uniform di [x0, x1]
y1 = runif(1, 0, mf) # Kandidat uniform di [0, mf]
y2 = f(u) # Nilai f(u)
if (y1 <= y2) { # Jika kandidat diterima
k = k + 1
hasil[k] = u # Simpan hasil
}
}
return(hasil) # Kembalikan hanya bilangan acak yang diterima
}
# Fungsi distribusi f(y)
f_y = function(y) {
ifelse(y <= 1, y, ifelse(y <= 1.5, 1, 0))
}
# Plot 1: Inverse Transform Method
hist(bilangan_acak, prob = TRUE, breaks = 20, col = "magenta", border = "white",
main = "(a) Inverse Transform Method", xlab = "y", ylab = "Density")
curve(ifelse(x <= 1, x, 1), add = TRUE, col = "red", lwd = 2)
# Menghasilkan 100 bilangan acak menggunakan Acceptance-Rejection Method
set.seed(123)
bil.acak = ARM_simplified(0, 1.5, 100, f_y)
bil.acak## [1] 1.4107009 0.8271525 1.4352500 1.0163560 0.3691316 1.3343090 0.9835587
## [8] 0.8160990 0.4337396 1.4445363 1.0360579 1.1376893 0.4772715 0.6205865
## [15] 0.2286671 1.1983873 0.8414220 1.3425680 0.9976728 0.5759545 1.2219601
## [22] 1.2150965 1.1915135 1.1317127 1.0652736 0.7129749 0.5276969 1.4774355
## [29] 1.3297036 0.4805599 1.1734415 0.7001686 0.8999834 1.3716573 0.6160347
## [36] 1.4029497 1.0808944 0.8782250 0.9718402 0.4615800 0.2313035 0.9288847
## [43] 1.0094986 0.7817036 1.2327082 1.4697329 0.4675533 1.2640940 0.3586499
## [50] 1.2711797 0.5818635 0.8579030 0.6671520 0.7534493 0.9749777 1.1105015
## [57] 0.6191192 0.9449596 1.2954662 1.0024270 0.5583571 1.3120235 1.2596516
## [64] 1.0624355 0.8915148 1.3697823 0.4112499 1.4784613 1.4059711 0.9020486
## [71] 0.5461378 0.2559679 0.7230639 0.5278330 1.2314270 1.0925916 0.7167681
## [78] 1.0473924 0.9275268 0.8131206 0.8202392 1.4510976 0.3858251 0.8895685
## [85] 0.7966056 1.3885619 1.0112803 0.7746673 0.5044968 1.3065651 1.1555011
## [92] 1.4578135 1.1378898 0.5948769 0.4466127 1.1313711 1.4732106 1.1768629
## [99] 1.1685988 0.9451978# Plot 2: Acceptance-Rejection Method
hist(bil.acak, prob = TRUE, breaks = 20, col = "pink", border = "white",
main = "(b) Acceptance Rejection Method", xlab = "y", ylab = "Density")
curve(ifelse(x <= 1, x, 1), add = TRUE, col = "red", lwd = 2)
# Plot kedua histogram secara bersamaan
par(mfrow = c(1, 2))
# Plot 1: Inverse Transform Method
hist(bilangan_acak, prob = TRUE, breaks = 20, col = "magenta", border = "white",
main = "(a) Inverse Transform Method", xlab = "y", ylab = "Density")
curve(ifelse(x <= 1, x, 1), add = TRUE, col = "red", lwd = 2)
# Plot 2: Acceptance-Rejection Method
hist(bil.acak, prob = TRUE, breaks = 20, col = "pink", border = "white",
main = "(b) Acceptance Rejection Method", xlab = "y", ylab = "Density")
curve(ifelse(x <= 1, x, 1), add = TRUE, col = "red", lwd = 2)