UAS PEMODELAN STATISTIKA DAN SIMULASI

1. Gunakan simulasi untuk menghitung nilai pendekatan dari E[X^2] untuk variabel random berdistribusi uniform jika X~U(0,1). Gunakan 10000 angka random

set.seed(123)
n <- 10000
x <- runif(n, min = 0, max = 1)

ekspektasi_x2 <- mean(x^2)
cat("Nilai ekspektasi simulasi:", ekspektasi_x2, "\n")
## Nilai ekspektasi simulasi: 0.3297405
ekspektasi_teoritis <- 1 / 3
cat("Nilai ekspektasi teoritis:", ekspektasi_teoritis, "\n")
## Nilai ekspektasi teoritis: 0.3333333

2. Buat masing-masing untuk n1 = 100, n2 = 1000, dan n3 = 10000 angka random dari N(100,15^2), lalu hitung rata-rata, simpangan baku, dan plot histogram. Berikan penjelasan hasil yang saudara peroleh.

set.seed(123)

analisis_simulasi <- function(n) {
  data <- rnorm(n, mean = 100, sd = 15)
  rata2 <- mean(data)
  simp_baku <- sd(data)
  
  # Hasil Statistik
  cat("Ukuran sampel =", n, "\n")
  cat("Rata-rata     =", round(rata2, 4), "\n")
  cat("Simpangan baku =", round(simp_baku, 4), "\n\n")
  
  # Plot histogram
  hist(data, breaks = 30, main = paste("Histogram N(100, 15^2), n =", n),
       xlab = "Nilai", col = "skyblue", border = "white")
}

# Simulasi n1
analisis_simulasi(100)
## Ukuran sampel = 100 
## Rata-rata     = 101.3561 
## Simpangan baku = 13.6922

# Simulasi n2
analisis_simulasi(1000)
## Ukuran sampel = 1000 
## Rata-rata     = 100.2862 
## Simpangan baku = 15.0689

# Simulasi n3
analisis_simulasi(10000)
## Ukuran sampel = 10000 
## Rata-rata     = 99.9634 
## Simpangan baku = 14.9831

Untuk ketiga sampel (n = 100, 1000, 10000), nilai rata-rata mendekati 100 dan simpangan baku mendekati 15, sesuai dengan parameter distribusi normal \(N(100, 15^2)\).

Semakin besar ukuran sampel, bentuk histogram semakin simetris dan menyerupai kurva normal sempurna.

Pada n = 100, histogram bisa tampak agak kasar, namun saat n meningkat, bentuk distribusi makin halus dan stabil.

3. Simulasikan distribusi binomial Bin(30, 0.25) sebanyak 5000 kali. Tentukan nilai estimasi P(X>=15)

set.seed(123)
x <- rbinom(5000, size = 30, prob = 0.25)

# Estimasi P(X >= 15)
estimasi_prob <- mean(x >= 15)
cat("Nilai Estimasi P(X >= 15):", estimasi_prob, "\n")
## Nilai Estimasi P(X >= 15): 0.0024

4. Hitunglah nilai s dengan membangkitkan data berdistribusi normal sebanyak data 1000, rata-rata 89, dan standar deviasi 10.

set.seed(123)

n <- 1000
mu <- 89
sigma <- 10

x <- rnorm(n, mean = mu, sd = sigma)

sum_x <- sum(x)
sum_x2 <- sum(x^2)
s <- sqrt((sum_x2 - (sum_x^2)/n) / (n - 1))
cat("Nilai s dari hasil simulasi adalah", s ,"\n" )
## Nilai s dari hasil simulasi adalah 9.91695

5. Sebuah bank memiliki 2 teller dengan kondisi berikut :

  • Nasabah datang rata-rata 5 orang per jam (dist Poisson)
  • Waktu pelayanan per nasabah rata-rata 8 menit (dist Eksponensial)
  • Bank buka 6 jam per hari

Kerjakan tugas berikut:

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.3.3
set.seed(123)

a. Buatlah simulasi antrial selama 20 hari kerja

set.seed(123)
lambda <- 5/60   
mu <- 1/8         
jam_kerja <- 360
n_hari <- 20


sim_bank <- function(n_teller) {
  total_tunggu <- 0
  total_nasabah <- 0
  total_sibuk <- 0
  antrian_harian <- matrix(0, n_hari, 6)
  
  for(hari in 1:n_hari) {
    
    n <- rpois(1, lambda * jam_kerja)
    kedatangan <- sort(runif(n, 0, jam_kerja))
    pelayanan <- rexp(n, mu)
    
   
    teller_selesai <- rep(0, n_teller)
    tunggu <- numeric(n)
    
    for(i in 1:n) {
      teller_idx <- which.min(teller_selesai)
      mulai <- max(kedatangan[i], teller_selesai[teller_idx])
      tunggu[i] <- mulai - kedatangan[i]
      teller_selesai[teller_idx] <- mulai + pelayanan[i]
    }
    
    
    total_tunggu <- total_tunggu + sum(tunggu)
    total_nasabah <- total_nasabah + n
    total_sibuk <- total_sibuk + sum(pmin(teller_selesai, jam_kerja))
 
    
    for(jam in 1:6) {
      start_t <- (jam-1) * 60
      end_t <- jam * 60
      dalam_jam <- kedatangan >= start_t & kedatangan < end_t
      if(sum(dalam_jam) > 0) {
        antrian_harian[hari, jam] <- mean(tunggu[dalam_jam] > 0)
      }
    }
  }
  
  list(
    rata_tunggu = total_tunggu / total_nasabah,
    utilisasi = total_sibuk / (n_teller * jam_kerja * n_hari),
    antrian_jam = colMeans(antrian_harian)
  )
}
hasil_2 <- sim_bank(2)
cat("Simulasi 20 hari dengan 2 teller selesai\n")
## Simulasi 20 hari dengan 2 teller selesai
cat("Total nasabah dilayani rata-rata:", round(rpois(1, lambda * jam_kerja)), "orang/hari\n")
## Total nasabah dilayani rata-rata: 32 orang/hari

b. Hitung berapa lama rata-rata nasabah menunggu

cat("Rata-rata waktu tunggu:", round(hasil_2$rata_tunggu, 2), "menit\n")
## Rata-rata waktu tunggu: 1.14 menit

c. Hitung berapa persen waktu teller sibuk bekerja

cat("Utilisasi teller:", round(hasil_2$utilisasi * 100, 1), "%\n")
## Utilisasi teller: 95.7 %

d. Bandingkan jika bank menambah 1 teller lagi - apakah efisien?

set.seed(123)
hasil_3 <- sim_bank(3)

perbandingan <- data.frame(
  Metrik = c("Waktu tunggu (menit)", "Utilisasi (%)", "Efisiensi"),
  Teller_2 = c(round(hasil_2$rata_tunggu, 2), 
               round(hasil_2$utilisasi * 100, 1),
               round(hasil_2$rata_tunggu / hasil_2$utilisasi, 2)),
  Teller_3 = c(round(hasil_3$rata_tunggu, 2),
               round(hasil_3$utilisasi * 100, 1), 
               round(hasil_3$rata_tunggu / hasil_3$utilisasi, 2))
)

print(perbandingan)
##                 Metrik Teller_2 Teller_3
## 1 Waktu tunggu (menit)     1.14     0.13
## 2        Utilisasi (%)    95.70    94.20
## 3            Efisiensi     1.19     0.14
penurunan <- (hasil_2$rata_tunggu - hasil_3$rata_tunggu) / hasil_2$rata_tunggu * 100
cat("\nPenurunan waktu tunggu:", round(penurunan, 1), "%\n")
## 
## Penurunan waktu tunggu: 88.7 %
cat("Kesimpulan:", ifelse(penurunan > 25, "EFISIEN - tambah teller", "KURANG EFISIEN"))
## Kesimpulan: EFISIEN - tambah teller

e. Buat grafik yang menunjukkan panjang antrian sepanjang hari

jam_data <- data.frame(
  Jam = 1:6,
  Teller_2 = hasil_2$antrian_jam,
  Teller_3 = hasil_3$antrian_jam
)

library(reshape2)
## Warning: package 'reshape2' was built under R version 4.3.3
jam_long <- melt(jam_data, id.vars = "Jam", variable.name = "Skenario", value.name = "Panjang")

ggplot(jam_long, aes(x = Jam, y = Panjang, fill = Skenario)) +
  geom_bar(stat = "identity", position = "dodge", alpha = 0.7) +
  labs(title = "Rata-rata Panjang Antrian per Jam",
       x = "Jam ke-", y = "Proporsi Waktu Ada Antrian") +
  theme_minimal() +
    scale_x_continuous(breaks = 1:6)

6. Bu Sari memiliki toko kelontong yang menjual beras dengan kondisi sebagai berikut :

  • Penjualan beras per hari : 8-15 karung (dist Uniform)
  • Stok awal : 100 karung
  • Setiap 5 hari sekali, Bu Sari pesan ulang 50 karung beras
  • Jika stok habis, kehilangan untung Rp50.000 per karung

Kerjakanlah tugas berikut :

a. Simulasikan stok beras selama 60 hari

library(ggplot2)
set.seed(123)
sim <- function(pesan=50, frek=5) {
  stok <- 100
  habis <- 0
  rugi <- 0
  data_stok <- c()
  data_jual <- c()
  
  for(i in 1:60) {
    jual <- round(runif(1, 8, 15))
    data_jual[i] <- jual
    
    if(stok < jual) {
      habis <- habis + 1
      rugi <- rugi + (jual - stok) * 50000
      stok <- 0
    } else {
      stok <- stok - jual
    }
    
    data_stok[i] <- stok
    
    if(i %% frek == 0) stok <- stok + pesan
  }
  
  list(stok=data_stok, habis=habis, rugi=rugi, jual=data_jual)
}

base <- sim()
cat("SIMULASI 60 HARI:")
## SIMULASI 60 HARI:
cat("\nRata-rata penjualan:", round(mean(base$jual), 1), "karung/hari")
## 
## Rata-rata penjualan: 11.6 karung/hari
cat("\nTotal terjual:", sum(base$jual), "karung")
## 
## Total terjual: 694 karung
cat("\nStok akhir:", tail(base$stok, 1), "karung")
## 
## Stok akhir: 0 karung
cat("\nHari kehabisan:", base$habis, "hari")
## 
## Hari kehabisan: 9 hari
cat("\nKerugian: Rp", base$rugi)
## 
## Kerugian: Rp 2200000

b. Hitung berapa kali toko kehabisan stok

cat("Kehabisan stok:", base$habis, "dari 60 hari")
## Kehabisan stok: 9 dari 60 hari
cat("\nPersentase:", round(base$habis/60*100, 1), "%")
## 
## Persentase: 15 %
cat("\nStok minimum:", min(base$stok))
## 
## Stok minimum: 0

c. Coba berbagai strategi: pesan 40, 50, atau 60 karung. Manakah yang terbaik?

hasil <- data.frame()
for(p in c(40, 50, 60)) {
  x <- sim(pesan=p)
  hasil <- rbind(hasil, c(p, x$habis, x$rugi))
}
names(hasil) <- c("Pesan", "Habis", "Rugi")
print(hasil)
##   Pesan Habis    Rugi
## 1    40    18 7900000
## 2    50     5 1850000
## 3    60     0       0
terbaik <- hasil$Pesan[which.min(hasil$Rugi)]
cat("\nTerbaik: pesan", terbaik, "karung")
## 
## Terbaik: pesan 60 karung

d. Coba ubah frekuensi pesan: setiap 3, 5, atau 7 hari. Manakah yang optimal?

hasil2 <- data.frame()
for(f in c(3, 5, 7)) {
  x <- sim(frek=f)
  hasil2 <- rbind(hasil2, c(f, x$habis, x$rugi))
}
names(hasil2) <- c("Frek", "Habis", "Rugi")
print(hasil2)
##   Frek Habis    Rugi
## 1    3     0       0
## 2    5     5 1950000
## 3    7    19 9050000
optimal <- hasil2$Frek[which.min(hasil2$Rugi)]
cat("\nOptimal: setiap", optimal, "hari")
## 
## Optimal: setiap 3 hari

e. Buat grafik pergerakan stok dan rekomendasi untuk bu Sari

# Simulasi optimal
opt <- sim(pesan=terbaik, frek=optimal)

# Plot
plot(1:60, base$stok, type="l", col="red", 
     xlab="Hari", ylab="Stok", main="Pergerakan Stok")
lines(1:60, opt$stok, col="blue")
legend("topright", c("Awal", "Optimal"), col=c("red", "blue"), lty=1)

# Rekomendasi
cat("\nREKOMENDASI:")
## 
## REKOMENDASI:
cat("\n- Pesan", terbaik, "karung setiap", optimal, "hari")
## 
## - Pesan 60 karung setiap 3 hari
cat("\n- Penghematan: Rp", base$rugi - opt$rugi)
## 
## - Penghematan: Rp 2200000
cat("\n- Kurangi kehabisan dari", base$habis, "ke", opt$habis, "hari")
## 
## - Kurangi kehabisan dari 9 ke 0 hari