membangkitkan bilangan acak sebanyak 10 dari distribusi normal menggunakan fungsi rnorm.

rnorm(10)

##  [1] -1.70163456 -2.85962485  0.85895446  0.70568803 -0.04435041  1.18994670
##  [7] -1.33597718  0.30369700  0.03896375 -1.38247749

Dengan mengatur set.seed sehingga menghasilkan bilangan acak yang diperoleh akan sama persis pada setiap pembangkitan.

set.seed(1)
rnorm(10)

##  [1] -0.6264538  0.1836433 -0.8356286  1.5952808  0.3295078 -0.8204684
##  [7]  0.4874291  0.7383247  0.5757814 -0.3053884

Membangkitkan bilangan acak normal dengan mean=5,dan σ 2 = 4

rnorm(10, mean=5, sd=2)

##  [1] 8.0235623 5.7796865 3.7575188 0.5706002 7.2498618 4.9101328 4.9676195
##  [8] 6.8876724 6.6424424 6.1878026

Membangkitkan bilangan acak dari distribusi uniform dengan fungsi runif.

runif(10)

##  [1] 0.8209463 0.6470602 0.7829328 0.5530363 0.5297196 0.7893562 0.0233312
##  [8] 0.4772301 0.7323137 0.6927316

Membangkitkan bilang acak dari distribusi uniform dengan nilai minimum 1 dan nilai maksimum 50

runif(10, 1, 100)

##  [1] 48.284343 86.259738 44.371614 25.234930  7.997226 10.847150 32.310899
##  [8] 52.344792 66.538503 41.276189

Mencari nilai peluang kumulatif peubah acak menggunakan distribusi normal.

p1 <- pnorm(1.645)
p1

## [1] 0.9500151

Mencari nilai kuantil dari peluang peubah acak

p1 <- pnorm(1.645)
p1

## [1] 0.9500151

Mencari nilai density peubah acak

a <- seq(-4, 4, length = 50)
da <- dnorm(a)
plot(a,da)

# Inverse-Transform method

Bangkitkan X~ Normal (10,1) sebanyak 100

r <- runif(100)
x <- qnorm(r, mean = 10, sd = 1) #invers dari fungsi sebaran kumulatif
hist(x, prob = TRUE, col="pink") # untuk melihat pola hasil datanya

sbx <- seq(7, 14, 0.01) # karena datanya dari 8 - 13, peningkatan garis by 0.01 supaya
smooth

## function (x, kind = c("3RS3R", "3RSS", "3RSR", "3R", "3", "S"), 
##     twiceit = FALSE, endrule = c("Tukey", "copy"), do.ends = FALSE) 
## {
##     if (!is.numeric(x)) 
##         stop("attempt to smooth non-numeric values")
##     if (anyNA(x)) 
##         stop("attempt to smooth NA values")
##     endrule <- match.arg(endrule)
##     rules <- c("copy", "Tukey")
##     if (is.na(iend <- pmatch(endrule, rules))) 
##         stop("invalid 'endrule' argument")
##     kind <- match.arg(kind)
##     if (startsWith(kind, "3RS") && !do.ends) 
##         iend <- -iend
##     else if (kind == "S") 
##         iend <- as.logical(do.ends)
##     type <- match(kind, c("3RS3R", "3RSS", "3RSR", "3R", "3", 
##         "S"))
##     smo <- .Call(C_Rsm, as.double(x), type, iend)
##     if (twiceit) {
##         r <- smooth(x - smo$y, kind = kind, twiceit = FALSE, 
##             endrule = endrule, do.ends = do.ends)
##         smo$y <- smo$y + r
##         if (!is.null(smo$iter)) 
##             smo$iter <- smo$iter + attr(r, "iter")
##         if (!is.null(smo$changed)) 
##             smo$changed <- smo$changed || attr(r, "changed")
##     }
##     if (is.ts(x)) 
##         smo$y <- ts(smo$y, start = start(x), frequency = frequency(x))
##     structure(smo$y, kind = kind, twiced = twiceit, iter = smo$iter, 
##         changed = smo$changed, endrule = if (startsWith(kind, 
##             "3")) 
##             rules[iend], call = match.call(), class = c("tukeysmooth", 
##             if (is.ts(x)) "ts"))
## }
## <bytecode: 0x0000029a9e47a048>
## <environment: namespace:stats>

lines(sbx, dnorm(sbx, mean = 10, sd = 1)) # gunakan denisty

Asumsikan Y adalah peubah acak berdistribusi eksponensial dengan λ = 2. Ingat kembali bahwa probability density function (pdf)-nya yaitu

p(y) = 2e−2y, untuk y > 0

Lakukan pembangkitan bilangan acak dengan metode inverse transformation.

n <- 1000
U <- runif(n)
X <- -log(1-U)/2
# plot
hist(X, freq=F, xlab='X',col="magenta", main='Exponential Random Variable')
curve(dexp(x, rate=2) , 0, 3, lwd=2, xlab = "", ylab = "", add = T)

Bangkitkan peubah acak yang berdistribusi Bernoulli dengan p = 0,4

n = 1000
p = 0.4
u = runif(n)
x = as.integer (u > 0.6)
barplot (table(x) , col = "Thistle")

Acceptance-rejection method

Bangkitkan peubah acak X~Beta (2,2) yang memiliki fungsi PDF f(x) = 6x(1 − x), 0 < x < 1 Dengan menggunakan Acceptance-Rejection Method sebanyak 1000 amatan. Kemudian gambarkan histogramnya.

Tetapkan peubah acak Y sebagai target dari sebaran beserta PDFnya g(y) Dalam kasus ini kita pilih Y~U(0,1) karena memiliki domain yang sama dengan peubah acak X yaitu 0 < x < 1. PDF dari Y yaitu g(y) = 1, 0 ≤ y ≤ 1
Bangkitkan peubah acak Y tersebut

#banyaknya amatan
n=1000
set.seed(10)
y=runif(n,0,1)

Bangkitkan u berdasarkan distribusi Uniform(0,1)

set.seed(30)
u=runif(n,0,1)

Hitung nilai c dengan mencari nilai maksimum dari f(t) g(t)

h<- function(y) 6*y*(1-y)
derivH <- optimize(f=h, interval= c(0,1), maximum=T)
derivH$objective

## [1] 1.5

Jika u ≤f(y)cg(y), jadikan Y sebagai X

f <- function(y) 6*y*(1-y)
g <- function(y) 1
nilaic <- derivH$objective
# kriteria penerimaan
kriteria <- u < f(y)/(nilaic*g(y))
head(kriteria)

## [1] TRUE TRUE TRUE TRUE TRUE TRUE

#jadikan Y sebagai X
x <- y[kriteria]
head(x)

## [1] 0.50747820 0.30676851 0.42690767 0.69310208 0.08513597 0.22543662

#banyaknya amatan
length(x)

## [1] 668

hist(x,main="beta(2,2)", col="pink")

#Direct Transformation

Contoh sebaran peubah acak khi-kuadrat berasal dari sebaran z dikuadratkan, sehingga untuk membangkitkan bilangan acak khi-kuadrat dengan db 1 dapat menggunakan z dikuadratkan.

set.seed(1)
y <- rnorm(1000)
x <- y^2
hist(x, prob = T)
sbx <- seq(0, 13, 0.01)
lines(sbx, dchisq(sbx, df = 1), col = "violet")

Transformasi dari distribusi uniform ke ditribusi normal dapat dilakukan dengan Transformasi Box-Muller.

n=1000
u1<-runif(n,0,1)
u2<-runif(n,0,1)
z1<-sqrt(-2*log(u1))*cos(2*pi*u2)
z2<-sqrt(-2*log(u1))*sin(2*pi*u2)
z<-c(z1,z2)
plot(density(z),xlab="x",ylab="p")
x<-seq(min(z),max(z),0.1)
y<-dnorm(x,0,1)
lines(x,y)

Tugas

Simulasikan peubah acak diskrit x berikut, dengan distribusinya sesuai dengan tabel berikut dengan metode Inverse-Transform method dengan n=1000

discrete_inverse_sampling <- function( prob ) { #definisi fungsi
  U  = runif(1)
  if(U <= prob[1]){
    return(1)
  }
  for(state in 2:length(prob)) {
    if(sum(prob[1:(state-1)]) < U && U <= sum(prob[1:state]) ) {
      return(state)
    }
  }
}

num_samples = 1000
prob = c(0.4, 0.2, 0.1, 0.1, 0.2)
names(prob) = c("0","2","3","7","10")
samples = runif(num_samples)
for(i in seq_len(num_samples) ) {
  samples[i] = discrete_inverse_sampling(prob)
}
sim_prob = table(samples) / sum(table(samples))
names(sim_prob) =  c("0","2","3","7","10")
barplot(prob, main='Fungsi Massa Probabilitas Sejati')

Disini, Barplot FMP sejati menunjukkan bahwa “0” memiliki nilai tertinggi yaitu sebesar 4.0.

barplot(sim_prob, main='Fungsi Massa Peluang Empiris')

Disini, FMP empiris menunjukkan bahwa “0” memiliki nilai tertinggi yaitu lebih dari 0.30.

Analisis : Setelah menghasilkan seribu angka acak, yang diambil dari distribusi Xi(0,2,3,7,10) dengan probabilitas P(X=xi) = (0.4, 0.2, 0.1, 0.1, 0.2), kita memperoleh sampel berikut: (0, 2, 3, 7, 10, 409, 194, 99, 101, 197). Kemudian, kami memvisualisasikan hasil ini dalam bentuk barplot untuk menganalisis fungsi masa probabilitas dan fungsi massa peluang empiris. Hasil yang ditemukan tidak berbeda jauh dengan distribusi awal, namun ada perbedaan kecil, terutama pada Xi(7) dalam fungsi massa peluang empiris yang memiliki nilai kurang dari 0.1.

Simulasikan peubah acak diskrit x berikut, dengan distribusinya sesuai dengan tabel berikut dengan metode Inverse-Transform method dengan n=1000

discrete_inverse_sampling <- function( prob ) {
  U  = runif(1)
  if(U <= prob[1]){
    return(1)
  }
  for(state in 2:length(prob)) {
    if(sum(prob[1:(state-1)]) < U && U <= sum(prob[1:state]) ) {
      return(state)
    }
  }
}

num_samples = 1000
prob = c(0.1, 0.2, 0.2, 0.2, 0.3)
names(prob) = c("0","1","2","3","4")
samples = runif(num_samples)
for(i in seq_len(num_samples) ) {
  samples[i] = discrete_inverse_sampling(prob)
}
sim_prob = table(samples) / sum(table(samples))
names(sim_prob) =  c("0","2","2","2","4")
barplot(prob, main='Fungsi Massa Probabilitas Sejati')

Disini, Barplot FMP sejati menunjukkan bahwa “4” memiliki nilai tertinggi yaitu sebesar 0.30.

barplot(sim_prob, main='Fungsi Massa Peluang Empiris')

Disini, Barplot FMP empiris menunjukkan bahwa “4” memiliki nilai tertinggi yaitu lebih dari 0.20.

Analisis : Hasil dari FMP sejati maupun empiris menunjukkan bahwa “4” memiliki nilai tertinggi. Perbedaannya adalah pada besaran sumbu y dimana pada FMP sejati berada pada posisi 0.30, sedangkan pada FMP empiris berada pada posisi lebih dari 0.20.

Bangkitkan bilangan acak dengan f(x) = 4e^−4x, denganx ∈ R. Gunakan metode Inverse- Transform method dengan n=1000

# Generate 1000 sampel dari distribusi
n <- 1000

samples <- numeric(n)
for (i in 1:n) {
  u <- runif(1)  # Generate variabel acak uniform antara 0 dan 1
  x <- (256/9) * u^(3/8)  # Menggunakan CDF invers
  samples[i] <- x
}

# Hasil sampel
hist(samples, breaks = 20, prob = TRUE, main = "Histogram of Samples", xlab = "X")

Analisis : Setelah menghasilkan seribu observasi dari variabel acak, dengan menggunakan metode invers transform dan menggunakan fungsi probabilitas f(x) = 4e^(-4x), di mana x adalah angka real, kami memperoleh hasil yang direpresentasikan dalam bentuk histogram. Dalam hasil observasi ini, terdapat tujuh nilai sampel, yaitu 0, 5, 10, 15, 20, 25, dan 30. Pada hasil ini, kita dapat melihat bahwa rentang nilai tertinggi terletak antara 25 dan 30, dengan tingkat kepadatan (density) sebesar 0.08.

Bangkitkan bilangan acak dengan f(x) =3/32x^5, dengan 0 < x < 2. Gunakan metode Inverse-Transform method dengan n=1000

# Define fungsi target PDF f(x)
f <- function(x) {
  return((3/2) * x^3 + (11/8) * x^2 + (1/6) * x + 1/2)
}

# Define fungsi proposal PDF g(x)
g <- function(x) {
  return(1)
}

# Hitung nilai M
x_values <- seq(0, 1, length.out = 1000)
M <- max(f(x_values) / g(x_values))

# Generate sampel menggunakan Acceptance-Rejection
n <- 1000
samples <- numeric(n)
count <- 0

while (count < n) {
  x_acak <- runif(1)  # Generate x dari g(x)
  u <- runif(1)      # Generate u dari distribusi uniform [0, 1]
  if (u <= f(x_acak) / (M * g(x_acak))) {
    samples[count + 1] <- x_acak
    count <- count + 1
  }
}

# Hasil sampel
hist(samples, breaks = 20, prob = TRUE, main = "Histogram of Samples", xlab = "X")

Analisis : Sama seperti dalam soal sebelumnya (soal nomor 3), pada soal nomor 4 ini, kita akan menampilkan sampel-sampel dalam bentuk histogram. Kali ini, kita memiliki 6 bilangan acak dengan nilai-nilai berikut: 0.0, 0.2, 0.4, 0.6, 0.8, dan 1.0. Dalam sampel ini, nilai tertinggi terletak pada 1.0, dengan nilai densitas yang melebihi 1.5.

Bangkitkan bilangan acak dengan metode Acceptance – rejection dengan fungsi PDF nya f(x) =3/2 x^3 + 11/8 x^2 + 1/6 x +1/2, 0 ≤ x ≤ 1

curve((3/2) * x^3 + (11/8) * x^2 + (1/6) * x + (1/12), 0, 1)

n <- 1000
j <- k <- 0
y <- numeric(n)
while (k < n) {
    u <- runif(1)
    j <- j + 1
    x <- runif(1)
    if (u < (((3/2) * (x^3)) + ((11/8) * (x^2)) + ((1/6) * x) + (1/12))/3.5 * (dunif(x, 
        0, 1))) {
        k <- k + 1
        y[k] <- x
    }
}
hist(y, prob = T)
sbx <- seq(0, 1, 0.01)
lines(sbx, ((3/2) * (sbx^3)) + ((11/8) * (sbx^2)) + ((1/6) * sbx) + (1/12))

Analisis : Metode Acceptance-Rejection, yang melibatkan beberapa langkah dalam proses pembuatan nilai acak sebagai berikut:

Pertama, kita tentukan peubah acak Y sebagai target dari distribusinya beserta fungsi kepadatan probabilitas (PDF) yang menggambarkan distribusi ini, yaitu g(y).
Selanjutnya, kita secara acak menghasilkan nilai-nilai dari peubah acak Y sesuai dengan distribusinya.
Selama proses ini, kita juga menghasilkan nilai u dari distribusi Uniform (0,1), yang sering dilambangkan sebagai U(0,1).
Kemudian, kita perhitungkan nilai c dengan mencari nilai maksimum dari hasil perkalian antara f(t) dan g(t), di mana t adalah peubah acak yang dihasilkan pada langkah (b).
Jika nilai u kurang dari atau sama dengan hasil perkalian antara f(y) dan c dengan g(y), maka kita akan menganggap Y sebagai nilai acak yang kita cari (X).
Langkah (e) diulang seiring dengan langkah (b), dan proses ini terus berlanjut, menciptakan sampel-sampel nilai acak yang akan mendekati distribusi probabilitas dari fungsi target yang telah ditentukan. Selama proses ini, kita menggambarkan kurva yang merepresentasikan fungsi target untuk keperluan visualisasi. Semakin banyak sampel yang dihasilkan (dengan n yang lebih besar), semakin baik perkiraan distribusi probabilitasnya.

Komstat 2

Jasmine Georgina Sekartaji_121450159

2023-10-04

Acceptance-rejection method

Tugas