UAS PSS

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

# Estimasi E[X²] untuk X ~ U(0,1)
set.seed(123)
x <- runif(10000, 0, 1)
E_X2 <- mean(x^2)
E_X2

## [1] 0.3297405

2. Buat masing-masing untuk𝑛1 = 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.

# Simulasi dari N(100, 15²)
set.seed(123)

simulasi_normal <- function(n) {
  data <- rnorm(n, mean = 100, sd = 15)
  mean_data <- mean(data)
  sd_data <- sd(data)
  
  hist(data, breaks = 30, col = "lightblue",
       main = paste("Histogram N =", n),
       xlab = "Nilai", ylab = "Frekuensi")
  
  cat("n =", n, "\n")
  cat("Rata-rata:", mean_data, "\n")
  cat("Simpangan baku:", sd_data, "\n\n")
}

simulasi_normal(100)

## n = 100 
## Rata-rata: 101.3561 
## Simpangan baku: 13.69224

simulasi_normal(1000)

## n = 1000 
## Rata-rata: 100.2862 
## Simpangan baku: 15.06891

simulasi_normal(10000)

## n = 10000 
## Rata-rata: 99.96336 
## Simpangan baku: 14.98314

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

# Simulasi Binomial Bin(30, 0.25)
set.seed(123)
sim_binom <- rbinom(5000, 30, 0.25)
prob_X_lebih_15 <- mean(sim_binom >= 15)
prob_X_lebih_15

## [1] 0.0024

4. Buatlah function untuk rumus berikut: sqrt((sum(x^2) - (sum(x)^2)/n) / (n - 1)) Hitung nilai s dengan membangkitkan data berdistribusi normal sebanyak data 1000, rata-rata 89, dan standar deviasi 10.

# Function untuk rumus simpangan baku
s_custom <- function(x) {
  n <- length(x)
  sqrt((sum(x^2) - (sum(x)^2)/n) / (n - 1))
}

set.seed(123)
x <- rnorm(1000, mean=89, sd=10)
s_custom(x)

## [1] 9.91695

5. (Simulasi Antrian di Bank) Sebuah bank memiliki 2 teller dengan kondisi berikut:

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

Kerjakanlah tugas berikut: a. Buatlah simulasi antrian selama 20 hari kerja b. Hitung berapa lama rata-rata nasabah menunggu c. Hitung berapa persen waktu teller sibuk bekerja d. Bandingkan jika bank menambah 1 teller lagi – apakah lebih efisien? e. Buat grafik yang menunjukkan panjang antrian sepanjang hari

5a. Buatlah Simulasi Antrian Selama 20 Hari Kerja dan 5b. Hitung Rata-rata Waktu Tunggu Nasabah

# Simulasi Antrian di Bank (Poisson & Eksponensial) Kondisi:
# λ = 5/jam, pelayanan: rata-rata 8 menit (μ = 7.5/jam), 2 teller
# 6 jam/hari × 20 hari kerja = 120 jam simulasi

# Simulasi antrian bank 20 hari kerja
library(simmer)

## Warning: package 'simmer' was built under R version 4.4.3

# Buat simulasi environment
env <- simmer("Bank Simulasi")

# Buat trajectory nasabah
traj <- trajectory("traj_nasabah") %>%
  seize("teller", 1) %>%
  timeout(function() rexp(1, rate = 1/8)) %>%   # Pelayanan rata-rata 8 menit
  release("teller", 1)

# Tambahkan resource dan generator sesuai kedatangan 5 orang per jam (setiap 12 menit)
env %>%
  add_resource("teller", capacity = 1) %>%
  add_generator("nasabah", traj, function() rexp(1, rate = 1/12)) %>%
  run(until = 7200)  # 6 jam × 20 hari = 7200 menit

## simmer environment: Bank Simulasi | now: 7200 | next: 7206.76351956559
## { Monitor: in memory }
## { Resource: teller | monitored: TRUE | server status: 0(1) | queue status: 0(Inf) }
## { Source: nasabah | monitored: 1 | n_generated: 611 }

# Ambil data kedatangan
data_arrival <- get_mon_arrivals(env, per_resource = TRUE)

# Hitung waktu tunggu: (end_time - activity_time) - start_time
data_arrival$waiting_time <- (data_arrival$end_time - data_arrival$activity_time) - data_arrival$start_time

# Tampilkan sebagian hasil
head(data_arrival[, c("start_time", "end_time", "activity_time", "waiting_time")])

##   start_time  end_time activity_time  waiting_time
## 1   19.96423  32.60213    12.6379000 -3.552714e-15
## 2   38.47329  44.32720     5.8539147  7.105427e-15
## 3   38.81971  56.90004    12.5728370  5.507492e+00
## 4  101.62154 103.35913     1.7375952  0.000000e+00
## 5  103.98794 112.01159     8.0236485  0.000000e+00
## 6  105.57822 112.62975     0.6181586  6.433375e+00

# Hitung rata-rata waktu tunggu
rata2 <- mean(data_arrival$waiting_time, na.rm = TRUE)
cat("5b. Rata-rata waktu tunggu nasabah:", round(rata2, 2), "menit\n")

## 5b. Rata-rata waktu tunggu nasabah: 14.69 menit

5c. Hitung Berapa Persen Waktu Teller Sibuk

# Ambil data resource
resource_data <- get_mon_resources(env)

# Urutkan berdasarkan waktu
resource_data <- resource_data[order(resource_data$time), ]

# Hitung durasi antar waktu
resource_data$duration <- c(diff(resource_data$time), 0)

# Hitung waktu sibuk total (server * durasi)
total_busy_time <- sum(resource_data$server * resource_data$duration)

# Total waktu simulasi
total_time <- 20000  # sesuai durasi run

# Hitung persentase waktu sibuk
persen_sibuk <- (total_busy_time / total_time) * 100

cat("✅ Persentase waktu teller sibuk:", round(persen_sibuk, 2), "%\n")

## ✅ Persentase waktu teller sibuk: 24.24 %

5d. Bandingkan jika bank menambah 1 teller lagi – apakah lebih efisien?

# Buat ulang environment baru
env_2teller <- simmer("Bank Simulasi - 2 Teller")

# Buat ulang trajectory yang sama
traj_2 <- trajectory("traj_nasabah") %>%
  seize("teller", 1) %>%
  timeout(function() rexp(1, rate = 1/8)) %>%
  release("teller", 1)

# Simulasi bank dengan 2 teller
env_2teller %>%
  add_resource("teller", capacity = 2) %>%
  add_generator("nasabah", traj_2, function() rexp(1, rate = 1/15)) %>%
  run(until = 20000)

## simmer environment: Bank Simulasi - 2 Teller | now: 20000 | next: 20003.2729956919
## { Monitor: in memory }
## { Resource: teller | monitored: TRUE | server status: 0(2) | queue status: 0(Inf) }
## { Source: nasabah | monitored: 1 | n_generated: 1280 }

# Ambil data kedatangan dan waktu tunggu
arrivals_2 <- get_mon_arrivals(env_2teller, per_resource = TRUE)
arrivals_2$waiting_time <- (arrivals_2$end_time - arrivals_2$activity_time) - arrivals_2$start_time

# Hitung rata-rata waktu tunggu
rata2_2teller <- mean(arrivals_2$waiting_time, na.rm = TRUE)
cat("✅ Rata-rata waktu tunggu dengan 2 teller:", round(rata2_2teller, 2), "menit\n")

## ✅ Rata-rata waktu tunggu dengan 2 teller: 0.39 menit

# Ambil data resource untuk 2 teller
res_2 <- get_mon_resources(env_2teller)
res_2 <- res_2[order(res_2$time), ]
res_2$duration <- c(diff(res_2$time), 0)

# Hitung total waktu sibuk
total_busy_2 <- sum(res_2$server * res_2$duration)

# Hitung persentase sibuk
persen_sibuk_2 <- (total_busy_2 / 20000) * 100
cat("✅ Persentase waktu teller sibuk (2 teller):", round(persen_sibuk_2, 2), "%\n")

## ✅ Persentase waktu teller sibuk (2 teller): 50.64 %

# Simulasi dengan 1 teller
env1 <- simmer("1 Teller")
traj1 <- trajectory("traj_nasabah") %>%
  seize("teller", 1) %>%
  timeout(function() rexp(1, 1/8)) %>%
  release("teller", 1)

env1 %>%
  add_resource("teller", capacity = 1) %>%
  add_generator("nasabah", traj1, function() rexp(1, 1/15)) %>%
  run(until = 20000)

## simmer environment: 1 Teller | now: 20000 | next: 20006.7350338103
## { Monitor: in memory }
## { Resource: teller | monitored: TRUE | server status: 0(1) | queue status: 0(Inf) }
## { Source: nasabah | monitored: 1 | n_generated: 1332 }

arr1 <- get_mon_arrivals(env1, per_resource = TRUE)
arr1$waiting_time <- (arr1$end_time - arr1$activity_time) - arr1$start_time
res1 <- get_mon_resources(env1)
res1 <- res1[order(res1$time), ]
res1$duration <- c(diff(res1$time), 0)
busy1 <- sum(res1$server * res1$duration)

# Simulasi dengan 2 teller
env2 <- simmer("2 Teller")
traj2 <- trajectory("traj_nasabah") %>%
  seize("teller", 1) %>%
  timeout(function() rexp(1, 1/8)) %>%
  release("teller", 1)

env2 %>%
  add_resource("teller", capacity = 2) %>%
  add_generator("nasabah", traj2, function() rexp(1, 1/15)) %>%
  run(until = 20000)

## simmer environment: 2 Teller | now: 20000 | next: 20040.895940186
## { Monitor: in memory }
## { Resource: teller | monitored: TRUE | server status: 0(2) | queue status: 0(Inf) }
## { Source: nasabah | monitored: 1 | n_generated: 1313 }

arr2 <- get_mon_arrivals(env2, per_resource = TRUE)
arr2$waiting_time <- (arr2$end_time - arr2$activity_time) - arr2$start_time
res2 <- get_mon_resources(env2)
res2 <- res2[order(res2$time), ]
res2$duration <- c(diff(res2$time), 0)
busy2 <- sum(res2$server * res2$duration)

# Buat tabel ringkasan
hasil <- data.frame(
  Skenario = c("1 Teller", "2 Teller"),
  Rata2_Waktu_Tunggu = c(mean(arr1$waiting_time, na.rm = TRUE),
                         mean(arr2$waiting_time, na.rm = TRUE)),
  Persen_Teller_Sibuk = c(busy1/20000 * 100, busy2/20000 * 100)
)

# Tampilkan hasil
print(hasil)

##   Skenario Rata2_Waktu_Tunggu Persen_Teller_Sibuk
## 1 1 Teller          8.8986927            51.90629
## 2 2 Teller          0.4765797            52.41769

5e. Buat Grafik Panjang Antrian Sepanjang Hari

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.4.2

# Ambil data dari simulasi terakhir (misal: env2 / 2 teller)
res_plot <- get_mon_resources(env2)

# Filter data antrian saja
res_queue <- subset(res_plot, resource == "teller")

# Buat plot panjang antrian
ggplot(res_queue, aes(x = time, y = queue)) +
  geom_line(color = "steelblue", size = 1) +
  labs(title = "📊 Panjang Antrian Sepanjang Hari",
       subtitle = "Simulasi Antrian Bank (2 Teller)",
       x = "Waktu (menit)",
       y = "Jumlah Antrian") +
  theme_minimal()

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

6. (Simulasi Toko kelontong) Bu Sari memliki toko kelontong yang menjual beras dengan kondisi berikut

Penjualan beras per hari: 8-15 karung (distribusi unifrom)
Stok awal: 100 karung
Bu Sari memesan 50 karung beras setiap 5 hari.
Jika stok habis, maka terjadi kerugian Rp50.000 per karung.
Simulasi dilakukan untuk 60 hari.

Tugas: a. Simulasikan stok 60 hari b. Hitung berapa kali kehabisan stok c. Coba strategi pemesanan 40, 50, 60 karung d. Coba ubah frekuensi pesan: tiap 3, 5, 7 hari e. Buat grafik pergerakan stok + rekomendasi

6a. Simulasikan stok beras selama 60 hari

set.seed(123)

# Parameter awal
hari <- 60
stok_awal <- 100
frekuensi_pesan <- 5
jumlah_pesan <- 50
kerugian_per_karung <- 50000

# Permintaan harian (uniform 8–15)
permintaan_harian <- sample(8:15, hari, replace = TRUE)

# Vektor hasil simulasi
stok <- numeric(hari)
stok[1] <- stok_awal - permintaan_harian[1]

for (i in 2:hari) {
  # Tambah stok jika hari pemesanan
  if ((i - 1) %% frekuensi_pesan == 0) {
    stok[i] <- stok[i-1] + jumlah_pesan
  } else {
    stok[i] <- stok[i-1]
  }

  # Kurangi stok sesuai permintaan
  stok[i] <- stok[i] - permintaan_harian[i]

  # Jika negatif, diset 0 untuk real stok
  stok[i] <- max(0, stok[i])
}

# Gabungkan jadi data frame
simulasi_stok <- data.frame(
  Hari = 1:hari,
  Permintaan = permintaan_harian,
  Stok = stok
)

head(simulasi_stok, 10)

##    Hari Permintaan Stok
## 1     1         14   86
## 2     2         14   72
## 3     3         10   62
## 4     4         13   49
## 5     5         10   39
## 6     6          9   80
## 7     7          9   71
## 8     8         13   58
## 9     9         10   48
## 10   10         12   36

6b. Hitung berapa kali kehabisan stok

jumlah_kehabisan <- sum(stok == 0)
cat("6b. Toko kehabisan stok sebanyak", jumlah_kehabisan, "hari\n")

## 6b. Toko kehabisan stok sebanyak 3 hari

6c. Coba berbagai strategi: 40, 50, 60 karung

simulasi_strategi <- function(jumlah_pesan) {
  stok <- numeric(hari)
  stok[1] <- stok_awal - permintaan_harian[1]
  for (i in 2:hari) {
    if ((i - 1) %% frekuensi_pesan == 0) {
      stok[i] <- stok[i-1] + jumlah_pesan
    } else {
      stok[i] <- stok[i-1]
    }
    stok[i] <- stok[i] - permintaan_harian[i]
    stok[i] <- max(0, stok[i])
  }
  sum(stok == 0)
}

strategi_hasil <- data.frame(
  Pesan = c(40, 50, 60),
  Hari_Habis = sapply(c(40, 50, 60), simulasi_strategi)
)

strategi_hasil

##   Pesan Hari_Habis
## 1    40         15
## 2    50          3
## 3    60          0

6d. Coba ubah frekuensi pesan: tiap 3, 5, 7 hari

simulasi_frekuensi <- function(frekuensi_pesan) {
  stok <- numeric(hari)
  stok[1] <- stok_awal - permintaan_harian[1]
  for (i in 2:hari) {
    if ((i - 1) %% frekuensi_pesan == 0) {
      stok[i] <- stok[i-1] + jumlah_pesan
    } else {
      stok[i] <- stok[i-1]
    }
    stok[i] <- stok[i] - permintaan_harian[i]
    stok[i] <- max(0, stok[i])
  }
  sum(stok == 0)
}

frekuensi_hasil <- data.frame(
  Frekuensi_Hari = c(3, 5, 7),
  Hari_Habis = sapply(c(3, 5, 7), simulasi_frekuensi)
)

frekuensi_hasil

##   Frekuensi_Hari Hari_Habis
## 1              3          0
## 2              5          3
## 3              7         19

6e. Grafik pergerakan stok + rekomendasi

library(ggplot2)

ggplot(simulasi_stok, aes(x = Hari, y = Stok)) +
  geom_line(color = "darkgreen", size = 1.2) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(title = "Grafik Pergerakan Stok Beras Bu Sari Selama 60 Hari",
       x = "Hari", y = "Stok (karung)") +
  theme_minimal()