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.

X <- runif(10000, min = 0, max = 1)

X_squared <- X^2

E_X_squared <- mean(X_squared)

cat("Pendekatan E[X^2] dengan simulasi:", E_X_squared, "\n")
## Pendekatan E[X^2] dengan simulasi: 0.3303049
set.seed(123)

simulasi_normal <- function(n, mu = 100, sigma = 15) {
    data <- rnorm(n, mean = mu, sd = sigma)
    rata_rata <- mean(data)
    simpangan_baku <- sd(data)

    hist(data, breaks = 30, main = paste("Histogram N(", mu, ",", sigma, ") n=", n), xlab = "Nilai", col = "skyblue", border = "white")
    abline(v = mu, col = "red", lwd = 2, lty = 2)
    abline(v = rata_rata, col = "blue", lwd = 2, lty = 2)

    cat("n = ", n, "\n")
    cat("Rata-rata = ", rata_rata, "\n")
    cat("Simpangan Baku Sampel = ", simpangan_baku, "\n")

}
simulasi_normal(n = 100)

## n =  100 
## Rata-rata =  101.3561 
## Simpangan Baku Sampel =  13.69224
simulasi_normal(n = 1000)

## n =  1000 
## Rata-rata =  100.2862 
## Simpangan Baku Sampel =  15.06891
simulasi_normal(n = 10000)

## n =  10000 
## Rata-rata =  99.96336 
## Simpangan Baku Sampel =  14.98314

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

set.seed(123)

simulasi <- rbinom(5000, size = 30, prob = 0.25)
proporsi_15 <- mean(simulasi >= 15)

cat("Estimasi P(X >= 15) dari simulasi adalah:", proporsi_15, "\n")
## Estimasi P(X >= 15) dari simulasi adalah: 0.0024

4. Buatlah function untuk rumus berikut: S = SQRT(Sigma x^2 - ((sigma x)^2/n) / n - 1) hitung nilai S dengan membangkitkan data berdistribusi normal sebanyak 1000, rata-rata = 89, dan standar deviasi 10

# Fungsi untuk menghitung S sesuai rumus
hitung_S <- function(x) {
  n <- length(x)
  S <- sqrt((sum(x^2) - (sum(x)^2)/n) / (n - 1))
  return(S)
}

# Bangkitkan data normal
set.seed(123)
data <- rnorm(1000, mean = 89, sd = 10)

# Hitung nilai S
S <- hitung_S(data)
cat("Nilai S hasil simulasi:", S, "\n")
## Nilai S hasil simulasi: 9.91695

5. (Simulasi Antrian Bank) sebuah bank memiliki 2 teller dengan kondisi berikut.

(a) Buatlah simulasi antrian selama 20 hari kerja

set.seed(123)

# Parameter 
lambda <- 5
jam_buka <- 6
n_teller <- 2
rata_pelayanan <- 8
hari_simulasi <- 20

# Fungsi simulasi satu hari
simulasi_hari <- function() {
    
    n_nasabah <- rpois(1, lambda * jam_buka)

    waktu_kedatangan <- sort(runif(n_nasabah, min = 0, max = jam_buka * 60))

    waktu_pelayanan <- rexp(n_nasabah, rate = 1/rata_pelayanan)

    selesai_teller <- rep(0, n_teller)
    waktu_mulai <- numeric(n_nasabah)
    waktu_selesai <- numeric(n_nasabah)
    waktu_tunggu <- numeric(n_nasabah)

    for (i in 1:n_nasabah) {
        teller_tersedia <- which.min(selesai_teller)
        waktu_mulai[i] <- max(waktu_kedatangan[i], selesai_teller[teller_tersedia])
        waktu_selesai[i] <- waktu_mulai[i] + waktu_pelayanan[i]
        waktu_tunggu[i] <- waktu_mulai[i] - waktu_kedatangan[i]
        selesai_teller[teller_tersedia] <- waktu_selesai[i]
    }

    return(list(
       n_nasabah = n_nasabah,
       waktu_tunggu = waktu_tunggu,
       waktu_pelayanan = waktu_pelayanan
    ))
}
set.seed(123)
hasil_simulasi <- lapply(1:hari_simulasi, function(x) simulasi_hari())
rata2_nasabah_perhari <- sapply(hasil_simulasi, function(x) x$n_nasabah)
rata2_waktu_tunggu_perhari <- sapply(hasil_simulasi, function(x) mean(x$waktu_tunggu))
cat("Rata-rata nasabah per hari:", mean(rata2_nasabah_perhari), "\n")
## Rata-rata nasabah per hari: 28.05
cat("Rata-rata waktu tunggu per hari:", mean(rata2_waktu_tunggu_perhari), "\n")
## Rata-rata waktu tunggu per hari: 1.009554

(b) Hitung berapa lama rata-rata nasabah menunggu

plot(rata2_nasabah_perhari, type = "b", pch = 19, col = "blue", lwd = 2, main = "Rata-rata Nasabah per Hari", xlab = "Hari", ylab = "Jumlah Nasabah")

plot(rata2_waktu_tunggu_perhari, type = "b", pch = 19, col = "red", lwd = 2, main = "Rata-rata Waktu Tunggu per Hari", xlab = "Hari", ylab = "Waktu Tunggu (menit)")

(c) Hitung berapa persen waktu teller sibuk bekerja

n_teller <- 2
waktu_kerja_perhari <- 6 * 60
total_waktu_kerja <- hari_simulasi * waktu_kerja_perhari

total_sibuk_teller <- rep(0, n_teller)

for (hari in 1:hari_simulasi) {
    n_nasabah <- hasil_simulasi[[hari]]$n_nasabah
    waktu_kedatangan <- sort(runif(n_nasabah, min = 0, max = waktu_kerja_perhari))
    waktu_pelayanan <- rexp(n_nasabah, rate = 1/rata_pelayanan)

    selesai_teller <- rep(0, n_teller)

    for (i in 1:n_nasabah) {
        teller_tersedia <- which.min(selesai_teller)
        waktu_mulai <- max(waktu_kedatangan[i], selesai_teller[teller_tersedia])
        waktu_selesai <- waktu_mulai + waktu_pelayanan[i]
        total_sibuk_teller[teller_tersedia] <- total_sibuk_teller[teller_tersedia] + waktu_pelayanan[i]
        selesai_teller[teller_tersedia] <- waktu_selesai
    }
}
# Hitung persentase waktu teller sibuk bekerja
persentase_sibuk <- total_sibuk_teller / total_waktu_kerja * 100

cat("Persentase waktu teller sibuk bekerja:", persentase_sibuk, "%\n")
## Persentase waktu teller sibuk bekerja: 31.15199 30.32007 %

(d) Bandingkan jika bank memiliki 1 teller lagi -apakah lebih efisien?

simulasi_hari_antrian <- function(n_teller) {
  n_nasabah <- rpois(1, lambda * jam_buka)
  waktu_kedatangan <- sort(runif(n_nasabah, min = 0, max = jam_buka * 60))
  waktu_pelayanan <- rexp(n_nasabah, rate = 1/rata_pelayanan)
  
  selesai_teller <- rep(0, n_teller)
  waktu_mulai <- numeric(n_nasabah)
  waktu_selesai <- numeric(n_nasabah)
  waktu_tunggu <- numeric(n_nasabah)
  panjang_antrian <- numeric(n_nasabah)
  
  for (i in 1:n_nasabah) {
    teller_tersedia <- which.min(selesai_teller)
    waktu_mulai[i] <- max(waktu_kedatangan[i], selesai_teller[teller_tersedia])
    waktu_selesai[i] <- waktu_mulai[i] + waktu_pelayanan[i]
    waktu_tunggu[i] <- waktu_mulai[i] - waktu_kedatangan[i]
    selesai_teller[teller_tersedia] <- waktu_selesai[i]
    
    # Hitung panjang antrian saat nasabah ke-i datang
    panjang_antrian[i] <- sum(selesai_teller > waktu_kedatangan[i])
  }
  
  return(list(
    n_nasabah = n_nasabah,
    waktu_tunggu = waktu_tunggu,
    waktu_pelayanan = waktu_pelayanan,
    panjang_antrian = panjang_antrian,
    waktu_kedatangan = waktu_kedatangan
  ))
}
set.seed(123)
hasil_2teller <- simulasi_hari_antrian(2)
hasil_3teller <- simulasi_hari_antrian(3)

cat("Rata-rata waktu tunggu 2 teller:", mean(hasil_2teller$waktu_tunggu), "menit\n")
## Rata-rata waktu tunggu 2 teller: 0.4839384 menit
cat("Rata-rata waktu tunggu 3 teller:", mean(hasil_3teller$waktu_tunggu), "menit\n")
## Rata-rata waktu tunggu 3 teller: 0 menit
cat("Rata-rata panjang antrian 2 teller:", mean(hasil_2teller$panjang_antrian), "\n")
## Rata-rata panjang antrian 2 teller: 1.538462
cat("Rata-rata panjang antrian 3 teller:", mean(hasil_3teller$panjang_antrian), "\n")
## Rata-rata panjang antrian 3 teller: 1.482759

(e)

plot(hasil_2teller$waktu_kedatangan, hasil_2teller$panjang_antrian, type = "s", col = "red", lwd = 2,
     xlab = "Menit sejak buka", ylab = "Panjang Antrian",
     main = "Panjang Antrian Sepanjang Hari (2 vs 3 Teller)")
lines(hasil_3teller$waktu_kedatangan, hasil_3teller$panjang_antrian, type = "s", col = "blue", lwd = 2)
legend("topright", legend = c("2 Teller", "3 Teller"), col = c("red", "blue"), lwd = 2)

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

(a) Simulasikan stok beras selama 60 hari

simulasi_stok_beras <- function(jumlah_pesan, interval_pesan, hari_simulasi = 60, stok_awal = 100, 
                                penjualan_min = 8, penjualan_max = 15, rugi_per_karung = 50000) {
  stok <- numeric(hari_simulasi + 1)
  stok[1] <- stok_awal
  rugi <- 0
  kehabisan <- 0
  stok_harian <- numeric(hari_simulasi)
  
  for (hari in 1:hari_simulasi) {
    # Penjualan hari ini
    penjualan <- sample(penjualan_min:penjualan_max, 1)
    # Jika stok cukup
    if (stok[hari] >= penjualan) {
      stok[hari + 1] <- stok[hari] - penjualan
    } else {
      # Jika stok tidak cukup, hitung rugi
      rugi <- rugi + (penjualan - stok[hari]) * rugi_per_karung
      kehabisan <- kehabisan + 1
      stok[hari + 1] <- 0
    }
    # Pesan ulang setiap interval_pesan hari
    if (hari %% interval_pesan == 0) {
      stok[hari + 1] <- stok[hari + 1] + jumlah_pesan
    }
    stok_harian[hari] <- stok[hari + 1]
  }
  return(list(stok_harian = stok_harian, rugi = rugi, kehabisan = kehabisan))
}
# Simulasi untuk 60 hari, pesan 50 karung setiap 5 hari
hasil_50 <- simulasi_stok_beras(50, 5)

# Tampilkan hasil simulasi stok harian
print(hasil_50$stok_harian)
##  [1] 86 74 60 46 87 74 61 47 32 68 58 49 37 24 59 44 29 17  3 50 39 28 15  4 50
## [26] 40 32 18  5 50 42 28 19  4 50 40 25 14  0 50 39 31 16  3 50 39 31 16  1 50
## [51] 40 25 14  0 50 38 27 16  5 50

(b) Hitung berapa kali toko kehabisan stok

# Hitung berapa kali kehabisan stok
cat("Selama 60 hari, toko kehabisan stok sebanyak", hasil_50$kehabisan, "kali.\n")
## Selama 60 hari, toko kehabisan stok sebanyak 10 kali.

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

set.seed(123)
hasil_40 <- simulasi_stok_beras(40, 5)
hasil_50 <- simulasi_stok_beras(50, 5)
hasil_60 <- simulasi_stok_beras(60, 5)
cat("Strategi pesan 40 karung:\n")
## Strategi pesan 40 karung:
cat("  Total kehabisan stok:", hasil_40$kehabisan, "hari\n")
##   Total kehabisan stok: 15 hari
cat("  Total rugi: Rp", hasil_40$rugi, "\n\n")
##   Total rugi: Rp 6400000
cat("Strategi pesan 50 karung:\n")
## Strategi pesan 50 karung:
cat("  Total kehabisan stok:", hasil_50$kehabisan, "hari\n")
##   Total kehabisan stok: 9 hari
cat("  Total rugi: Rp", hasil_50$rugi, "\n\n")
##   Total rugi: Rp 2850000
cat("Strategi pesan 60 karung:\n")
## Strategi pesan 60 karung:
cat("  Total kehabisan stok:", hasil_60$kehabisan, "hari\n")
##   Total kehabisan stok: 0 hari
cat("  Total rugi: Rp", hasil_60$rugi, "\n\n")
##   Total rugi: Rp 0

(d) Coba ubah frekuensi pesan: 40, 50, atau 60 karung. Manakah yang terbaik?

# Simulasi untuk 3, 5, dan 7 hari
hasil_3hari <- simulasi_stok_beras(jumlah_pesan = 60, interval_pesan = 3)
hasil_5hari <- simulasi_stok_beras(jumlah_pesan = 60, interval_pesan = 5)
hasil_7hari <- simulasi_stok_beras(jumlah_pesan = 60, interval_pesan = 7)
# Tampilkan hasil kehabisan stok
cat("Pesan setiap 3 hari: kehabisan stok", hasil_3hari$kehabisan, "kali\n")
## Pesan setiap 3 hari: kehabisan stok 0 kali
cat("Pesan setiap 5 hari: kehabisan stok", hasil_5hari$kehabisan, "kali\n")
## Pesan setiap 5 hari: kehabisan stok 0 kali
cat("Pesan setiap 7 hari: kehabisan stok", hasil_7hari$kehabisan, "kali\n")
## Pesan setiap 7 hari: kehabisan stok 10 kali

(E) Buat grafik pergerakan stok dan rekomendasi untuk Bu Sari

# simulasi pergerakan
plot(hasil_40$stok_harian, type = "l", col = "red", lwd = 2, ylim = c(0, max(hasil_60$stok_harian)),
     xlab = "Hari", ylab = "Stok Beras", main = "Pergerakan Stok Beras (60 Hari)")
lines(hasil_50$stok_harian, col = "blue", lwd = 2)
lines(hasil_60$stok_harian, col = "green", lwd = 2)
legend("topright", legend = c("Pesan 40", "Pesan 50", "Pesan 60"), col = c("red", "blue", "green"), lwd = 2)

# Simulasi frekuensi pesan
plot(hasil_3hari$stok_harian, type = "l", col = "red", lwd = 2, ylim = c(0, max(hasil_3hari$stok_harian, hasil_5hari$stok_harian, hasil_7hari$stok_harian)),
     xlab = "Hari", ylab = "Stok Beras", main = "Pergerakan Stok Beras (Frekuensi Pesan Berbeda)")
lines(hasil_5hari$stok_harian, col = "blue", lwd = 2)
lines(hasil_7hari$stok_harian, col = "green", lwd = 2)
legend("topright", legend = c("Pesan 3 hari", "Pesan 5 hari", "Pesan 7 hari"), col = c("red", "blue", "green"), lwd = 2)